## Abstract

Much is known about the nonlinear resonant response of mechanical systems, but methods for the systematic design of structures that optimize aspects of these responses have received little attention. Progress in this area is particularly important in the area of micro-systems, where nonlinear resonant behaviour is being used for a variety of applications in sensing and signal conditioning. In this work, we describe a computational method that provides a systematic means for manipulating and optimizing features of nonlinear resonant responses of mechanical structures that are described by a single vibrating mode, or by a pair of internally resonant modes. The approach combines techniques from nonlinear dynamics, computational mechanics and optimization, and it allows one to relate the geometric and material properties of structural elements to terms in the normal form for a given resonance condition, thereby providing a means for tailoring its nonlinear response. The method is applied to the fundamental nonlinear resonance of a clamped–clamped beam and to the coupled mode response of a frame structure, and the results show that one can modify essential normal form coefficients by an order of magnitude by relatively simple changes in the shape of these elements. We expect the proposed approach, and its extensions, to be useful for the design of systems used for fundamental studies of nonlinear behaviour as well as for the development of commercial devices that exploit nonlinear behaviour.

## 1. Introduction

The nonlinear resonant behaviour of structures is quite well understood for the case of single-mode resonance and for nonlinear interactions between a few modes. The phenomena associated with these resonances include amplitude-dependent natural frequencies, instabilities and bifurcations, coexistence of multiple steady-states, etc. There are abundant analytical and simulation studies on these topics, and plenty of experimental evidence to support the validity of simple models for describing these phenomena. In many applications, these behaviours are viewed as curiosities, generally to be avoided.

However, recent advances in the fields of nonlinear micro- and nano-resonators [1] have contributed to the development of new MEMS devices that use nonlinear dynamics, e.g. the hardening and softening behaviour of individual modes and energy transfer between modes, for applications, including energy harvesting [2], atomic force microscopy [3], mass detection [4], inertial sensing [5], passive frequency division [6] and frequency stabilization [7]. In most previous studies, devices have been designed that demonstrate the behaviour of interest, but more recently researchers have started to explore the possibility of tailoring systems for targeted responses. For single-mode resonators, tuning the nonlinear resonances with hardening and softening behaviour is considered by Cho *et al.* [8], where the orientation of a nano-tube connecting two cantilever structures was changed to obtain hardening or softening behaviour, and Dai *et al.* [9], where the impact of topology variation on the hardening behaviour was investigated in dynamic tests of an energy harvesting device. For coupled-mode resonators with internal resonances, Tripathi & Bajaj [10] presented a procedure to synthesize families of structures with two commensurable natural frequencies (such as 1:2 or 1:3 relations). By contrast, Pedersen [11] investigated structural optimization for minimizing the possibility of internal resonance by avoiding integer relations between the frequencies. To the authors’ knowledge, structural manipulation and optimization of the essential nonlinear effects that capture the energy transfer between coupled modes has not been reported. These effects are naturally described by the normal form associated with a given resonance condition [12], in which the coefficients of particular nonlinear terms dictate the behaviour at the given resonance; we refer to these as the *essential* nonlinearities of the system. The key to optimizing nonlinear resonant responses is in the manipulation of these coefficients by structural modifications.

Gradient-based structural optimization is a powerful tool in design optimization. In previous research, it has been applied in nonlinear structural dynamics to optimize the amplitude of transient responses [13,14], periodic responses [15,16] and also eigenfrequencies of deformed structures with geometrical nonlinearities [17]. For a clamped–clamped beam, optimal shape design associated with nonlinear frequency responses, including primary resonances and super-harmonic resonances, was studied by Dou & Jensen [18] using the incremental harmonic balance method. Another related work is the optimal shape design of comb fingers for electrostatic forces as investigated by Ye *et al.* [19], which produced combinations of linear, quadratic and cubic driving-force profiles. Shaped comb drives have been fabricated, tested and widely used in many applications, including tunable resonators [20] and large-displacement parametric resonators [21].

Here, a gradient-based optimization method is presented that allows one to tailor the nonlinearities associated with certain nonlinear resonances for a given resonator configuration. First, the coefficients associated with the *essential* nonlinearity for a given resonance are explicitly derived from nonlinear finite-element models; this involves finite-element discretization, explicit modal expansion and distillation of the normal form coefficients. Then, the quantity that captures the *essential* nonlinearity in the reduced-order model is manipulated by structural optimization. For example, the hardening and softening behaviour associated with the first flexural mode of a clamped–clamped beam is captured by the coefficient of a certain cubic term in a single-mode reduced-order model, while for internal resonances the main nonlinearity is a nonlinear modal coupling term that promotes energy exchange between modes.

For a micro-resonator with complex geometry, characterization based on nonlinear finite-element models offers a powerful tool to derive reduced-order models. This procedure has been investigated for the analysis of nonlinear vibration of piezoelectric layered beams with applications to NEMS [22]. A more general theory for finite-element models with geometric nonlinearity was provided by Touzé *et al.* [23] with applications to reduced-order models by using ‘Mixed Interpolation of Tensorial Components’ (MITC) shell elements. Here we adopt this method for characterization of nonlinear resonators, which is applicable at macro and micro scales. Our focus is on nonlinearities arising from finite deformation kinematics of structures, but the methodology can be adapted to other situations for which forces can be described in terms of a potential, such as electrostatic loads, and for which the material nonlinearity, i.e. a nonlinear stress–strain relation, is expressed in terms of polynomial functions.

The article is organized as follows. First, the characterization theory is introduced in §2, where the coefficients for nonlinear stiffness and nonlinear modal coupling, for both the Hamiltonian and the reduced-order differential equation model, are explicitly derived from the finite-element model. Continuing in §2, the sensitivity of the nonlinear modal coupling coefficients and the general optimization problem are formulated. Examples given in §§3 and 4 provide some discussion, conclusion and directions for future work.

## 2. Analysis and optimization

### (a) Characterization theory

In order to demonstrate the methodology, we consider systems for which the dominant nonlinear effects arise from non-inertial conservative effects and we neglect the nonlinear inertial terms that might be important in, for example, cantilevered structures. This includes planar frame structures for which mid-plane stretching is the primary nonlinearity. For convenience in implementation, the derivation is formulated in matrix-vector form instead of tensor form. The full set of coefficients for individual mode nonlinear stiffness and for nonlinear modal coupling are first derived for the Hamiltonian of the system, and then the attendant essential coefficients in the reduced-order model are obtained in a straightforward manner.

The characterization theory is derived for general finite-element models with quadratic and cubic nonlinearities arising from nonlinear strain–displacement relations. With the finite-element discretization of the continuous structure, the Hamiltonian for the system, i.e. the sum of kinetic energy *T* and potential energy *U*, can be expressed as
2.1where **M** is the mass matrix of the finite-element model, **u** is the global vector of nodal displacements, ** ϵ** and

**are element-wise vectors of strain and stress components, respectively, and**

*σ**V*

_{e}indicates that the volume integration is performed within the element

*e*. Assuming a linear strain–stress relation, the potential energy

*U*can be written as 2.2where is a symmetric constitutive matrix. We may now divide the strain into a linear and nonlinear part as 2.3where

**B**

_{0}is the linear strain–displacement matrix and

**B**

_{1}(

**u**

^{e}) is a function of the element-wise vector of nodal displacements

**u**

^{e}. Substituting equation (2.3) into equation (2.2) and using the symmetry of , the potential energy

*U*can be divided into three parts, representing its expansion at leading orders, as 2.4The first step in setting up the modal equations is to solve the corresponding linear eigen-problem 2.5where

**K**is the linear stiffness matrix of the finite-element model, to yield a set of modal frequencies and mode shapes (

*ω*

_{p},

*Φ*_{p}); these are normalized w.r.t.

**M**such that and . We then express the displacements as a superposition of

*N*

_{m}linear modes as 2.6where

*q*

_{p}are the modal coordinates and are the element-wise mode shape vectors extracted from the global vector. Substituting equation (2.6) into equation (2.4), we obtain the strain energy expressed in terms of mode shapes and modal coordinates as 2.7where the linear modal coupling coefficients and the nonlinear modal coupling coefficients, and , are explicitly expressed as 2.8Since we use the linear eigenmodes normalized with respect to the mass matrix, we have for

*i*≠

*j*and . With these modal coupling coefficients, we can now write the Hamiltonian for the system in equation (2.1) in modal coordinates (

*q*

_{i},

*p*

_{i}) (where ) as 2.9where the relation has been used. With the Hamiltonian of the system, it is straightforward to derive a set of ordinary differential equations in modal coordinates as 2.10where

*p*=1,…,

*N*

_{m},

*q*

_{p}is the modal coordinate corresponding to

*Φ*_{p}, and we have introduced the modal force

*f*

_{p}obtained from a projection of the load vector

**f**in the full finite-element model onto the

*p*th mode, i.e. , as well as modal damping ratios (damping as a ratio of critical damping) expressed as

*ξ*

_{p}. The modal coupling coefficients and are 2.11It is noted that the differential equation representation of equation (2.10) with the modal coupling coefficients in equation (2.11) is equivalent to the formulation obtained using the principle of virtual work in previous studies [23].

It is well known that a large number of linear modes may be required to accurately describe the behaviour of nonlinear resonators. On the other hand, the hardening/softening behaviour of a single-mode resonator is well predicted by a single second-order governing equation with its dynamics projected onto a single nonlinear normal mode. This advanced reduced-order model is achieved by using normal form theory, linked to nonlinear normal modes, as described in [24,25]. Using the theory of normal forms, our modal coordinates (*q*_{p},*p*_{p}) can be transformed into the curved, normal coordinates (*R*_{p},*S*_{p}) (where ) through a nonlinear transformation of coordinates, written in general form as
2.12For details, the reader is referred to [24,25]. For clarity, the reduced-order model and the corresponding nonlinear modal coupling coefficients in the curved coordinates of the nonlinear normal modes are further discussed in the examples for a single-mode resonator and a coupled-mode resonator with internal resonance.

In this work, we apply the framework to structures modelled by beam elements [18]. However, it should be noted that the characterization and optimization procedure in this paper works independent of the choice of elements.

### (b) Optimization and sensitivity analysis

In this section, we present the general formulation for optimizing an objective function that depends on nonlinear coefficients of interest for a given model. For example, we can maximize/minimize the hardening behaviour in a clamped–clamped beam by maximizing/minimizing the coefficient of the cubic nonlinearity.

For generality, we consider an objective function *c* that may be an explicit function of the nonlinear coefficients, as well as the eigenvectors and associated eigenvalues, and formulate our optimization problem as the following minimization problem:
2.13where the subscripts *i*, *j*, *k*, *l*, *p*=1,…,*N*_{m}. The optimization problem is subjected to a set of constraints associated with the beam shape parametrization: *h*_{e} is the thickness of a beam element which is bounded by in the optimization and *ρ*_{e} is the normalized design variable bounded between 0 and 1. In practice, the lower bound is dictated by fabrication tolerances, and the upper bound is used to keep the beam relatively slender. For coupled -mode resonators, such as the one treated in the second example of this paper, we impose an additional constraint such that the ratio of two associated eigenvalues stays in a sufficiently small neighbourhood of *n*_{1}*ω*_{p1}=*n*_{2}*ω*_{p2}, where, for internal resonance conditions, *n*_{1} and *n*_{2} are selected integers and *p*_{1} and *p*_{2} are the orders of the two modes of interest.

For efficient structural optimization, we will use a gradient-based approach. The sensitivity of the objective function can be calculated by using direct differentiation [26,27], but a more efficient approach for many design variables is the adjoint method [28] where we merely need to solve *N*_{m} groups of adjoint equations. To derive the adjoint equation, the objective function is first rewritten with adjoint variables **λ**_{p} and *η*_{p} as
2.14where it is noted that the terms in the two sets of parentheses that appear in the appended term both vanish identically. Differentiation of the objective function with respect to design variable *ρ*_{e} is then expressed as
2.15Since the direct computation of d*Φ*_{p}/d*ρ*_{e} and d*ω*_{p}/d*ρ*_{e} is computationally expensive, the adjoint variables are selected such that the coefficients of d*Φ*_{p}/d*ρ*_{e} and d*ω*_{p}/d*ρ*_{e} vanish. This leads to the adjoint equations in terms of **λ**_{p} and *η*_{p} as
2.16and
2.17where we use the symmetry of **M** and **K**, and
2.18In this case, the sensitivity of the objective function writes
2.19

It is noted that the two quantities and are assembled from the corresponding element-wise quantities and . For convenience of computational implementation, the adjoint equations are expressed in matrix form as 2.20Based on the objective function and calculated sensitivities, the update of design variables is performed using the mathematical optimization tool called the method of moving asymptotes (MMA) [29], which solves a series of convex approximating subproblems. The algorithm has proved efficient for large-scale structural optimization. A new system analysis is then performed with the updated design variables. These steps are repeated until the design variables no longer change within some prescribed small tolerance. To summarize the procedure, a flow chart of the proposed optimization is displayed in figure 1. We now demonstrate the approach with examples.

## 3. Examples

Two examples are presented. The first is to maximize/minimize the cubic nonlinearity of the fundamental mode in a nonlinear resonator consisting of a clamped–clamped beam, a common element in MEMS applications. The second example shows how one can maximize the essential modal coupling nonlinearity in a T-bar frame with a 2:1 internal resonance, a structure that has been proposed as an MEMS frequency divider.

### (a) Single-mode resonator

A crucial feature of a nonlinear resonator is the hardening and softening behaviour associated with a given vibration mode. For a lightly damped single-mode resonator, we will focus on the hardening and softening behaviour of its free responses, i.e. without damping and external loads. This problem has been investigated from a structural dynamics/finite-element perspective (e.g. [24]). For its applicability in general structures, including symmetric structures like clamped–clamped beams and asymmetric structures like arches or shells, the model development includes both quadratic and cubic terms. For single-mode resonator, the reduced-order model based on a single linear mode is written as
3.1with and . The frequency–amplitude relation is derived as
3.2where *a* is the amplitude of the corresponding linear mode, and the effective coefficient *Γ* is
3.3

A more accurate model for the hardening/softening behaviour of these structures can be created with the dynamics projected onto a single nonlinear normal mode and the corresponding equation is then written as
3.4with the new coefficients and given as explicit functions of nonlinear modal coupling coefficients and , [24]. Based on equation (3.4), the frequency–amplitude relationship is derived as
3.5where *A* is the amplitude of the corresponding nonlinear normal mode in curved, normal coordinates and the effective coefficient *Γ** is
3.6as approximated by averaging methods. Note that in this case *Γ** is not only linked to the *p*th linear mode explicitly through *ω*_{p} and and implicitly through which is in and , but also linked to other linear modes, whose contributions are taken into account through and .

We find that for the flexural mode of a clamped–clamped beam with geometric nonlinearity from mid-plane stretching, the dominating coefficient is and we can therefore simplify our optimization problem considerably by approximating *Γ** as
3.7since , due to the symmetry. Based on the coefficient *Γ** in equation (3.7) and omitting the constant factor, the objective function is now selected as
3.8

where the plus/minus sign corresponds to minimizing/maximizing the hardening behaviour, respectively. This simplification is possible since we consider the fundamental mode of a clamped–clamped beam, which has a strictly hardening nonlinearity, that is . Substituting the objective function *c* in equation (3.8) into equation (2.19), its sensitivity with respect to design variables *ρ*_{e} is obtained as
3.9with adjoint variables **λ**_{p} and *η*_{p} solved from adjoint equation in equation (2.20) with
3.10where
3.11The beam has a fixed length *L* of 300 μm and a fixed out-plane width of 20 μm. The initial design has a uniform in-plane thickness of 4 μm and is discretized with 400 beam elements as described in [18]. During shape optimization, the in-plane thickness *h* is varied to tailor the cubic nonlinearity in the reduced-order model. We set , and . The material properties are assumed for Si, that is mass density *ρ*=2329 kg m^{−3}, and Young’s modulus *E*=170 GPa. The vibration modes of the initial design and two optimized designs are shown in figures 2, 3 and 4, respectively. Evolution of the objective function and shapes obtained during the evolution are shown in figures 5 and 6. In these optimizations, the objective function is increased by a factor of 13 and reduced by a factor of 4, respectively. The optimized designs are in accordance with the results in [18], obtained using the incremental harmonic balance method, where we found the nonlinear strain energy due to mid-plane stretching reaches its local maximum around and , which is precisely where the optimized structures are altered most significantly relative to their general thickness. Furthermore, the eigenfrequency of the first flexural mode decreases during optimization of maximizing the cubic nonlinearity, and increases during optimization of minimizing the cubic nonlinearity. This follows from the fact that the structure is made generally thinner when maximizing the nonlinearity, so that the critical sections can be made relatively thick, and the opposite trend occurs when minimizing the nonlinearity. It should be emphasized that, no guarantee can be made that the obtained designs are global optima and the found solutions will in general depend on the chosen initial design. However, repeated optimization runs reveal that the formulation is quite robust and we believe only a marginal improvement in performance could be possibly found with another design.

### (b) Coupled-mode resonator with internal resonance

We present an example of two-to-one internal resonance for which the normal form has a single important inter-modal coupling term. This example is motivated by an MEMS frequency divider that makes use of internal resonance, for which a general theory is presented in [6]. The structure considered is also similar to other proposed MEMS devices developed for filtering [30,31]. In these structures, two vibrational modes can exchange energy during free vibration. For convenience, we will refer to modes *p*_{1} and *p*_{2} as modes 1 and 2, respectively, in our discussion. The divider device consists of two localized modes with *ω*_{p2}≈2*ω*_{p1}, which provides a division by two in frequency, as measured in mode 1, when mode 2 is driven near its resonance. The term ‘localized mode’ indicates that the dominant vibration associated with the mode occurs in a localized part of the structure, even though the entire structure is generally involved in the modal vibration. For instance, for the T-bar structure shown in figure 7, the vibration of mode 1 is localized in the vertical beam and the vibration of mode 2 is localized in the horizontal beam. In operation, a harmonic load with frequency close to *ω*_{2} is applied to drive mode 2 into resonance. When the amplitude of mode 2 is sufficiently large, it will induce the vibration of mode 1 due to the parametric pumping, wherein the transverse vibration of the horizontal beam provides an axial force in the vertical beam, which in turn induces transverse motion of the vertical beam when the frequency of the horizontal beam is approximately twice that of the vertical beam. The response of the structure is governed by a model in which these two modes are coupled through resonant terms. In application, *n* such elements are linked to achieve division by 2^{n}; division by eight has been experimentally achieved [32]. The reduced-order dynamic model with its conservative dynamics projected onto two nonlinear normal modes in curved, normal coordinates is written as
3.12and
3.13where the explicit expressions of and are given in [24]. The key term of interest is that associated with the *essential* modal coupling coefficients and , since they are the terms in the normal form that describes the resonant nonlinear coupling terms that promote energy exchange between the modes. These modal coupling coefficients can also be observed from the governing equation with its dynamics projected onto two linear modes, which is written as
3.14where the most important coefficients are and , between which there is an exact relation . The reason for this exact relation is that the coefficients of *q*_{1}*q*_{2} and in equation (3.14) arise from differentiation of the term in the Hamiltonian with respect to *q*_{1} and *q*_{2}, where . In fact, this example is particularly attractive for demonstrating the present approach since energy exchange between the modes can be described (to leading order) by this single nonlinear term. Also, as observed in the equation for *q*_{1} in the reduced-order model, *β* is essentially the amplitude of the parametric excitation provided to mode 1 from mode 2, given by 2*β*_{1}*q*_{2}*q*_{1}, and is therefore related to the stability region of the linear parametric resonance of *q*_{1}. The term in the equation for *q*_{2} captures the back-coupling of the driven mode onto the driving mode, as required for passive coupling, which is beneficial in frequency dividers [32]. The combined effect of this coupling is the possibility of energy transfer between the modes when there is a two-to-one internal resonance.

The optimization problem for this resonance is therefore formulated as
3.15where *ϵ*=0.001 and additional constraints are imposed as in equation (2.13). The sensitivity is computed with equation (2.19) written as
3.16where . It is noted that there are two groups of adjoint variables, **λ**_{1} and *η*_{1}, **λ**_{2} and *η*_{2}, corresponding to the two vibration modes. These two groups of adjoint variables are solved using adjoint equation in equation (2.20) with *p*=1, 2.

For a specific example, the length of the horizontal beam is taken to be 300 μm and the length of the vertical beam is taken to be 195.5 μm, so that *ω*_{2}≈2*ω*_{1}.The lengths of the two beams are fixed during the optimization. The initial in-plane thickness is uniformly 4 μm along both beams, and the in-plane thickness is bounded between 2 μm and 6 μm during the optimization. The material properties are the same as in example 1. An optimized design and its two important vibration modes are displayed in figure 8. Evolution of the objective function and optimized designs over iterations are displayed in figure 9.

In order to demonstrate the effects of tuning and optimizing the coupling nonlinearity in this example, we carry out simulations of the structure for different shapes encountered during the optimization iteration process. With this internal resonance and zero damping, free vibrations of the system will exhibit energy exchange between the modes in a beating type response, whose features depend crucially on the magnitude of the coupling coefficient and the initial conditions, as well as the other system parameters. For each structural shape and initial amplitude of the starting mode, the exchange of energy has a particular beat period; this period will be shorter for higher starting amplitudes and for larger values of the coupling coefficient, since both these effects enhance the nonlinear modal coupling. Figure 10*a* shows a typical response obtained from the finite-element model of the final (optimized) shape, obtained by initiating the response using the second linear mode at a moderate amplitude. Figure 10*b* shows the normalized beat period for several values of the initial energy for three different designs, that is for three different values of the coupling coefficient. The predicted trends are evident as both the coupling coefficient and initial amplitude are varied. Note that it would be also worthwhile to compare results from the finite-element model with simulations and analysis of the reduced two mode model, but this requires a more detailed study of the effects of all nonlinear coefficients that vary during the optimization process.

Evolution of the eigenfrequencies of linear vibration modes 1 and 2 encountered over iterations of the optimization process is displayed in figure 11*a*. Other measures of interest for this system are (i) the degree of spatial energy localization in the vertical beam of the first vibration mode (note that from symmetry the second vibration mode has perfect localization) and (ii) the effectiveness of the horizontal beam in parametrically pumping the vertical beam in the second mode. The localization of the first mode is measured by the maximum transverse amplitude of the horizontal beam divided by the maximum transverse amplitude of the vertical beam. Likewise, the pumping effectiveness of the horizontal beam in the second mode is measured by the ratio of the transverse vibration at the midspan of the horizontal beam to the maximum transverse vibration of the same beam, which occurs near the quarter spans. The results in figure 11*b* show that the localization ratio decreases during optimization, which indicates improved localization, and the pumping ratio increases, which indicates enhanced coupling of the two modes. The intuition of the optimized design is that the axial deformation along the vertical beam increases with a larger end mass and a thinner cross section, both of them occur in the optimized design.

## 4. Conclusion

A systematic procedure is proposed for characterization and optimal design of nonlinear resonators from finite-element models. The method makes use of a means of computing terms in normal forms from nonlinear finite-element models and uses sensitivities to update design parameters to vary a given objective function, subject to constraints. The proposed procedure is demonstrated on the optimization of the nonlinear modal stiffness of a single-mode resonator and the optimization of an essential inter-mode coupling term in a resonator with two internally resonant modes. It is shown that quite simple shape alterations can provide significant changes in these nonlinear effects. Thus, this approach offers an important tool for the development of structures with desired nonlinear behaviour, which has particular applicability in micro-systems. Of course, these nonlinear features have a significant effect on the forced response of the system, and this is where most applications will benefit from tailoring the nonlinear response of the system. Also, it is interesting to note that by minimizing the nonlinearity in a single-mode structure, one can achieve a larger range of vibration amplitudes over which the system behaves according to a linear model; this may be of significant benefit in applications where one wishes to maximize the linear dynamic range of the device, for example, in frequency sources (oscillators) in which resonating elements are employed in a feedback loop to generate a periodic signal.

Current efforts are focused on experimentally verifying the approach and applying it to specific devices for frequency generation and conversion. Also, although only planar structures with geometrically stiffness nonlinearities are considered in the paper, it is possible to extend the methods to include material nonlinearities, inertial nonlinearities, nonlinear electrostatic forces and three-dimensional physics; these are especially simple for effects that can be modelled by a potential. One can also allow for more variability in structural shapes, for example, by allowing beam lengths to change, and one can impose additional frequency constraints, which could be important in applications. In addition, varying the shape of the structure changes all coefficients in the reduced two mode model, and these in turn affect the response, so more specific applications may require the use of objective functions that depend on several nonlinear coefficients. Also, while shape optimization is considered here, the philosophy of the present approach can also be applied using topology optimization, which may offer more elaborate structural configurations with better results than can be achieved by shape optimization. Work along these lines for structural models based on hyper-elastic materials has recently been reported [33].

## Authors' contributions

All four authors contributed to the general overview in section 1 and the conclusions in section 4. S.D. and B.S.S. were primarily responsible for the analytical development in section 2 and the examples in section 3. S.D. prepared the first draft of the manuscript. S.W.S. and J.S.J. provided general supervision and guidance throughout the effort.

## Competing interests

We declare we have no competing interests.

## Funding

This project was supported by ERC Starting grant no. 279529 INNODYN, US NSF grant no. 1234067 and DARPA-MTO-DEFYS.

## Acknowledgements

Special thanks to O. Shoshani, K. Turner, D. Czaplewski, D. Lopez, T. Kenny and C. Touzé for inspiring applications and useful discussions. S.D. thanks Otto Mønsteds Fond and Thomas B. Thriges Fond for supporting his research stay in Michigan State University in 2014.

## Appendix A

The linear and nonlinear formulation of the beam element used here can be found in article [18]. Its nonlinear modal coupling coefficients in the Hamiltonian can be written as
and
where *E* is Young’s modulus, *A* is cross-sectional area, **u**^{e} (the longitudinal displacements) and **w**^{e} (the transverse displacements) are obtained from the element-wise eigenvector *Φ*^{e} projected onto the element coordinates, the subscripts *i*, *j* and *k* indicate the order of the eigenvector, and with *L* denoting the element length and **K**_{g} are expressed as

## Footnotes

One contribution of 11 to a theme issue ‘A field guide to nonlinearity in structural dynamics’.

- Accepted March 13, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.