Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch article

Structural optimization for nonlinear dynamic response

Suguang Dou

Suguang Dou

Department of Mechanical Engineering, Technical University of Denmark, Kgs. Lyngby, Denmark

Google Scholar

Find this author on PubMed

,
B. Scott Strachan

B. Scott Strachan

Department of Mechanical Engineering, Michigan State University, East Lansing, MI, USA

Department of Electrical and Computer Engineering, Michigan State University, East Lansing, MI, USA

Google Scholar

Find this author on PubMed

,
Steven W. Shaw

Steven W. Shaw

Department of Mechanical Engineering, Michigan State University, East Lansing, MI, USA

Department of Physics and Astronomy, Michigan State University, East Lansing, MI, USA

[email protected]

Google Scholar

Find this author on PubMed

and
Jakob S. Jensen

Jakob S. Jensen

Department of Mechanical Engineering, Technical University of Denmark, Kgs. Lyngby, Denmark

Department of Electrical Engineering, Technical University of Denmark, Kgs. Lyngby, Denmark

Google Scholar

Find this author on PubMed

Published:https://doi.org/10.1098/rsta.2014.0408

    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

    Display Formula
    2.1
    where 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 Ve 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
    Display Formula
    2.2
    where Inline Formula is a symmetric constitutive matrix. We may now divide the strain into a linear and nonlinear part as
    Display Formula
    2.3
    where B0 is the linear strain–displacement matrix and B1(ue) is a function of the element-wise vector of nodal displacements ue. Substituting equation (2.3) into equation (2.2) and using the symmetry of Inline Formula, the potential energy U can be divided into three parts, representing its expansion at leading orders, as
    Display Formula
    2.4
    The first step in setting up the modal equations is to solve the corresponding linear eigen-problem
    Display Formula
    2.5
    where 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 Inline Formula and Inline Formula. We then express the displacements as a superposition of Nm linear modes as
    Display Formula
    2.6
    where qp are the modal coordinates and Inline Formula 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
    Display Formula
    2.7
    where the linear modal coupling coefficients Inline Formula and the nonlinear modal coupling coefficients, Inline Formula and Inline Formula, are explicitly expressed as
    Display Formula
    2.8
    Since we use the linear eigenmodes normalized with respect to the mass matrix, we have Inline Formula for ij and Inline Formula. With these modal coupling coefficients, we can now write the Hamiltonian for the system in equation (2.1) in modal coordinates (qi,pi) (where Inline Formula) as
    Display Formula
    2.9
    where the relation Inline Formula has been used. With the Hamiltonian of the system, it is straightforward to derive a set of ordinary differential equations in modal coordinates as
    Display Formula
    2.10
    where p=1,…,Nm, qp is the modal coordinate corresponding to Φp, and we have introduced the modal force fp obtained from a projection of the load vector f in the full finite-element model onto the pth mode, i.e. Inline Formula, as well as modal damping ratios (damping as a ratio of critical damping) expressed as ξp. The modal coupling coefficients Inline Formula and Inline Formula are
    Display Formula
    2.11
    It 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 (qp,pp) can be transformed into the curved, normal coordinates (Rp,Sp) (where Inline Formula) through a nonlinear transformation of coordinates, written in general form as

    Display Formula
    2.12
    For 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:

    Display Formula
    2.13
    where the subscripts i, j, k, l, p=1,…,Nm. The optimization problem is subjected to a set of constraints associated with the beam shape parametrization: he is the thickness of a beam element which is bounded by Inline Formula in the optimization and ρe is the normalized design variable bounded between 0 and 1. In practice, the lower bound Inline Formula is dictated by fabrication tolerances, and the upper bound Inline Formula 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 n1ωp1=n2ωp2, where, for internal resonance conditions, n1 and n2 are selected integers and p1 and p2 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 Nm groups of adjoint equations. To derive the adjoint equation, the objective function is first rewritten with adjoint variables λp and ηp as

    Display Formula
    2.14
    where 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
    Display Formula
    2.15
    Since 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
    Display Formula
    2.16
    and
    Display Formula
    2.17
    where we use the symmetry of M and K, and
    Display Formula
    2.18
    In this case, the sensitivity of the objective function writes
    Display Formula
    2.19

    It is noted that the two quantities Inline Formula and Inline Formula are assembled from the corresponding element-wise quantities Inline Formula and Inline Formula. For convenience of computational implementation, the adjoint equations are expressed in matrix form as

    Display Formula
    2.20
    Based 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.
    Figure 1.

    Figure 1. A flow chart of the proposed structural optimization for nonlinear dynamic response.

    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

    Display Formula
    3.1
    with Inline Formula and Inline Formula. The frequency–amplitude relation is derived as
    Display Formula
    3.2
    where a is the amplitude of the corresponding linear mode, and the effective coefficient Γ is
    Display Formula
    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

    Display Formula
    3.4
    with the new coefficients Inline Formula and Inline Formula given as explicit functions of nonlinear modal coupling coefficients Inline Formula and Inline Formula, [24]. Based on equation (3.4), the frequency–amplitude relationship is derived as
    Display Formula
    3.5
    where A is the amplitude of the corresponding nonlinear normal mode in curved, normal coordinates and the effective coefficient Γ* is
    Display Formula
    3.6
    as approximated by averaging methods. Note that in this case Γ* is not only linked to the pth linear mode explicitly through ωp and Inline Formula and implicitly through Inline Formula which is in Inline Formula and Inline Formula, but also linked to other linear modes, whose contributions are taken into account through Inline Formula and Inline Formula.

    We find that for the flexural mode of a clamped–clamped beam with geometric nonlinearity from mid-plane stretching, the dominating coefficient is Inline Formula and we can therefore simplify our optimization problem considerably by approximating Γ* as

    Display Formula
    3.7
    since Inline Formula, due to the symmetry. Based on the coefficient Γ* in equation (3.7) and omitting the constant factor, the objective function is now selected as
    Display Formula
    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 Inline Formula. Substituting the objective function c in equation (3.8) into equation (2.19), its sensitivity with respect to design variables ρe is obtained as

    Display Formula
    3.9
    with adjoint variables λp and ηp solved from adjoint equation in equation (2.20) with
    Display Formula
    3.10
    where
    Display Formula
    3.11
    The 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 Inline Formula, and Inline Formula. 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 Inline Formula and Inline Formula, 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.
    Figure 2.

    Figure 2. Initial design of a clamped–clamped beam and its linear vibration mode. The colour in online version indicates the vibration amplitude. (Online version in colour.)

    Figure 3.

    Figure 3. Optimized design for maximizing the cubic nonlinearity of a clamped–clamped beam and its linear vibration mode. The colour in online version indicates the vibration amplitude. (Online version in colour.)

    Figure 4.

    Figure 4. Optimized design for minimizing the cubic nonlinearity of a clamped–clamped beam and its linear vibration mode. The colour in online version indicates the vibration amplitude. (Online version in colour.)

    Figure 5.

    Figure 5. Evolution of the objective function and shapes encountered during the optimization process. The vertical axis is the absolute value of the objective function divided by its initial value and the horizontal axis is the iteration number. (Online version in colour.)

    Figure 6.

    Figure 6. Evolution of the objective function and shapes encountered during the optimization process. The vertical axis is the absolute value of the objective function divided by its initial value and the horizontal axis is the iteration number. (Online version in colour.)

    (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 p1 and p2 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 2n; 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

    Display Formula
    3.12
    and
    Display Formula
    3.13
    where the explicit expressions of Inline Formula and Inline Formula are given in [24]. The key term of interest is that associated with the essential modal coupling coefficients Inline Formula and Inline Formula, 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
    Display Formula
    3.14
    where the most important coefficients are Inline Formula and Inline Formula, between which there is an exact relation Inline Formula. The reason for this exact relation is that the coefficients of q1q2 and Inline Formula in equation (3.14) arise from differentiation of the term Inline Formula in the Hamiltonian Inline Formula with respect to q1 and q2, where Inline Formula. 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 q1 in the reduced-order model, β is essentially the amplitude of the parametric excitation provided to mode 1 from mode 2, given by 2β1q2q1, and is therefore related to the stability region of the linear parametric resonance of q1. The Inline Formula term in the equation for q2 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.
    Figure 7.

    Figure 7. Initial design and the two coupled linear vibration modes obtained using COMSOL modal analysis. (a) Linear vibration mode 1 and (b) linear vibration mode 2 and ω2≈2ω1. The colour in online version indicates the vibration amplitude. (Online version in colour.)

    The optimization problem for this resonance is therefore formulated as

    Display Formula
    3.15
    where ϵ=0.001 and additional constraints are imposed as in equation (2.13). The sensitivity is computed with equation (2.19) written as
    Display Formula
    3.16
    where Inline Formula. 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.

    Figure 8.

    Figure 8. Optimized design for maximizing the absolute value of the essential modal coupling coefficient and the two coupled linear vibration modes obtained using COMSOL modal analysis. (a) Linear vibration mode 1 and (b) linear vibration mode 2 and ω2≈2ω1. The colour in online version indicates the vibration amplitude. (Online version in colour.)

    Figure 9.

    Figure 9. Evolution of the objective function and shapes encountered during the optimization process. The vertical axis is the absolute value of the objective function divided by its initial value and the horizontal axis is the iteration number. (Online version in colour.)

    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 10a 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 10b 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.

    Figure 10.

    Figure 10. (a) Typical time simulation with added envelope curves showing beats, for the final design with initial first mode amplitude of 7.5×10−11 (using eigenvectors normalized by the mass matrix). (b) Data points and fitted exponential curves for the beat period, normalized by the period of the first linear mode, versus the initial amplitude of the first mode, for designs corresponding to iteration numbers 0 (initial design, top curve), 30 (middle curve) and 250 (final design, bottom curve). (Online version in colour.)

    Evolution of the eigenfrequencies of linear vibration modes 1 and 2 encountered over iterations of the optimization process is displayed in figure 11a. 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 11b 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.

    Figure 11.

    Figure 11. (a) Evolution of eigenfrequencies of linear vibration mode 1 and 2 encountered during the optimization process. (b) Evolution of pumping ratio and localization ratio encountered during the optimization process. (Online version in colour.)

    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 Inline Formula can be written as

    Display Formula
    and
    Display Formula
    where E is Young’s modulus, A is cross-sectional area, ue (the longitudinal displacements) and we (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 Inline Formula and Kg are expressed as
    Display Formula

    Footnotes

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

    Published by the Royal Society. All rights reserved.

    References