## Abstract

This paper investigates the development of a novel framework and its implementation for the nonlinear tuning of nano/microresonators. Using geometrically exact mechanical formulations, a nonlinear model is obtained that governs the transverse and longitudinal dynamics of multilayer microbeams, and also takes into account rotary inertia effects. The partial differential equations of motion are discretized, according to the Galerkin method, after being reformulated into a mixed form. A zeroth-order shift as well as a hardening effect are observed in the frequency response of the beam. These results are confirmed by a higher order perturbation analysis using the method of multiple scales. An inverse problem is then proposed for the continuation of the critical amplitude at which the transition to nonlinear response characteristics occurs. Path-following techniques are employed to explore the dependence on the system parameters, as well as on the geometry of bilayer microbeams, of the magnitude of the dynamic range in nano/microresonators.

## 1. Introduction

The selection of design parameters for nano- and microelectromechanical system (N/MEMS) devices is challenged by high parameter sensitivity [1] and uncertainty [2–4] in both the manufacturing process and operating environment. This is further exacerbated by the relatively high costs for manufacturing and testing microdevices *in situ*. Such limitations call for physical models that can predict the microdevice behaviour, within certain ranges of operation similar to actual applications, as accurately as possible. However, when the behaviour of N/MEMS devices is nonlinear, either intrinsically or extrinsically, the selection of proper design parameters from the corresponding model is, at best, accomplished on a case-by-case basis [5].

Schemes for parameter selection are often based on optimization of measurable properties of device operation [6], and/or the identification of desirable operating conditions (e.g. pull-in/lift-off [7]). Examples of strategies based on optimization include minimization of electrical sensitivity in piezoresistive cantilever sensors [6], maximization of tip deflection in curved-electrode actuators [8] and minimization of force variations in the presence of uncertainties in manufacturing using changes to device geometry [9]. As an alternative to optimization, in many applications of electrostatic microdevices, parameters are chosen to correspond to operation near static [10–12] or dynamic [13] pull-in. Similarly, the design of some ultrasensitive microsensors relies on parameter choices corresponding to cyclic-fold/saddle-node bifurcation points in their frequency response [14,15].

Numerical continuation offers a powerful tool for the identification of special operating conditions, as well as the solution to optimization problems. For example, the common simulation approach for detecting a pull-in point, which is realized as a fold of the solution manifold with respect to parameter space, is to sweep along the manifold until no convergence (or very slow convergence) is achieved. By contrast, a numerical continuation algorithm may be formulated to locate the fold point as a regular solution point [7,16], as well as to continue past the fold point along a branch of unstable responses. Moreover, brute-force simulation schemes for solving the optimization problem in the design of N/MEMS devices may be very time-consuming [17]. Proper formulation of a continuation problem can be used to accelerate the search for optimal solutions in the desired parameter space [18].

In this paper, our focus is on studying axially restrained beams [19]. Such structures find widespread application in micro- and nanoelectromechanical devices, from sensing [20] to detection [21] and communication systems [22]. These applications rely on predictive models of high fidelity that can be employed in robust device design. For reliable design of micro/nanoscale devices, low measurement signal-to-noise ratios require models that take into account all the phenomenological aspects that govern the device dynamics. This includes accounting for the different intrinsic nonlinear signatures of micro- and nanoscale beams [23].

Nonlinearities of the hardening type occur in the context of axially restrained planar beams. The hardening behaviour is due to geometric nonlinearities associated with stretching of the beam’s axis. In this case, the beam tension couples with the curvature in a way that becomes pronounced at high amplitudes of oscillation. This effect has been studied extensively [24,25], mostly inspired by the work of Mettler [26]. A second type of nonlinearity has a softening effect and is mainly associated with nonlinear inertia forces and the mass distribution in the specific geometric design of the beam as well as its boundary conditions.

In addition to intrinsic sources of nonlinearities in beams, the particular physics relied upon by the design of an electromechanical device can also contribute to nonlinear behaviour. This is the case, for example, with electrostatic actuation which depends in a strongly nonlinear manner on the beam deflection. Both intrinsic and extrinsic nonlinearities in beams have been used, favourably, in designing sensor [27–29], switch [30] and amplifier [31] applications, in which the nonlinear effects play a crucial role (for a review, see [32]). Despite the successful approach of exploiting nonlinearities at the design stage of electromechanical devices, nonlinear behaviours are, in general, considered as undesirable. They add uncertainty (due to multistability and coexistence of attractors) to the prediction of the system performance, and great care is often taken to reduce their influence [33].

For low amplitudes of excitation, the response of an axially restrained beam is essentially that of a linear structure, with a unique oscillatory behaviour for each excitation frequency. For sufficiently large amplitudes, the hardening effect becomes pronounced and bistability occurs over a portion of the frequency response. The desire to operate within the linear regime imposes a natural upper limit on the input excitation amplitude and, consequently, on the response amplitude. In many practical situations, including the application to nanoresonators, low signal-to-noise ratios further impose a lower bound on the input excitation, in order to result in a response that is distinguishable above the noise floor. Designs that increase the difference (known as the *dynamic range*) between the maximum response amplitude, while operating in the linear regime, and the noise floor are highly desirable.

The idea of tuning the nonlinearities in microelectromechanical devices has been investigated in the literature [34–36]. More recently, nonlinear tuning has been used to increase the dynamic range of electrostatically driven resonators [37–39] (but see also [40] for an example of relying on modal coupling). The need for analytical or computational schemes for finding critical values at which the transition to the nonlinear regime occurs has been highlighted by several authors [33,38,41,42]. However, these works are based on highly simplified assumptions both in the nonlinear model of the microdevice and in the approximation of the critical values.

It is the goal of this study to explore the influence of the beam geometry for a bilayer microbeam on (i) the magnitude of intrinsic nonlinear effects, and (ii) the dynamic range. To this end, in §2, we develop a distributed-mass geometrically nonlinear model of the microbeam dynamics. In §3, we apply a perturbation method to obtain an approximation for the free response of a single-layer nonlinear microbeam, which serves as a theoretical benchmark in the validation of the numerical scheme used in the subsequent discussion. The computational analysis is based on a finite-element discretization of a mixed weak formulation of the partial differential equations that govern the beam dynamics, which is described in §4. This method, in the context of an inverse problem that captures the dynamic range of the microbeam as a constraint, is used in §5 to explore this range across a subset of values of several design parameters.

## 2. Model derivation

Consider a microbeam with rectangular cross section and constant width. In this study, we restrict attention to planar motions [43].

### (a) Kinematics

Let *r*^{0}(*s*):=*s**e*_{1} be the position vector describing the base line of a stress-free straight beam parametrized by the horizontal coordinate *s*. Denote by ** r**(

*s*,

*t*):=

*s*

*e*_{1}+

*u*(

*s*,

*t*)

*e*_{1}+

*v*(

*s*,

*t*)

*e*_{2}the deformed configuration of the base line. Furthermore, as shown in figure 1, introduce a rotated reference triad to represent the current orientation of the material cross section at

*s*, given by

*r*_{s}=(1+

*u*

_{s})

*e*_{1}+

*v*

_{s}

*e*_{2}, such that

*s*denotes partial differentiation with respect to

*s*. Here,

*ν*−1,

*η*, and

*μ*:=

*θ*

_{s}denote the

*beam elongation*,

*shear strain*and

*flexural curvature*, respectively. Finally, let

*angular velocity vector*.

### (b) Dynamics

The balance of the linear and angular momentum in local form now yields the equations of motion [43]
** n**:=

*N*

*b*_{1}+

*Q*

*b*_{2}includes the tension

*N*and the shear force

*Q*,

**:=**

*m**M*

*b*_{3}is the bending moment,

**is the external force per unit reference length,**

*f***:=**

*c**c*

*b*_{3}is the external couple per unit reference length,

*ϱA*is the mass of the beam per unit reference length,

*A*is the cross-sectional area,

*ϱ*

**is the cross-sectional second mass moment of inertia per unit reference length, and**

*J**ϱ*

**is the vector of**

*i**first mass moment of inertia*, i.e.

*ϱ*

**:=**

*i**ϱI*

*b*_{2}, where

*ϱI*is the cross-sectional first mass moment of inertia. The intrinsic component form of the equations of motion becomes

*f*

_{1}and

*f*

_{2}are the components of

**along the**

*f*

*b*_{1}and

*b*_{2}directions, respectively.

### (c) Bilayer microbeams

Consider the special case of a bilayer beam whose first layer is the beam substrate with non-zero, constant width and thickness throughout the beam span. Note that the horizontal axis in figure 1 passes through the centreline of this substrate, and that the corresponding first area moment of inertia, *I*_{1}, equals zero about *b*_{3}. Here and henceforth, subscripts denote the beam layer and not the axis about which the area moments are computed.

Let *E*_{1} and *E*_{2} denote the Young moduli for the substrate and the second layer, respectively. The equations of motion then take the non-dimensional form
*E*_{1}*J*_{1}/ℓ^{2}, and lengths by ℓ, where *ϱ*_{i}*A*_{i} and *ϱ*_{i}*J*_{i} are the mass per unit length and the mass moment of inertia, respectively, for the *i*th layer. Here,
*I*_{2}=−*A*_{2}*h*/2, where *h*=*h*_{1}+*h*_{2}, and *h*_{1} and *h*_{2} denote the thicknesses of the substrate and second layer, respectively.

Polysilicon microbeams are often fabricated through an etching process that removes an underlying sacrificial layer to create a structural element with clamped ends. In this paper, we model such a structure by appending the clamped–clamped boundary conditions
*α* and λ are the slenderness ratio and the coefficient for the contribution of rotary effects, respectively.

### (d) Constitutive relationships

To obtain constitutive relationships between *N* and *M* and the beam’s generalized strains, consider the strain–displacement relationship for the elongation *ϵ*(*z*)=*ν*−1+*z* *μ* and the linearly elastic constitutive law *σ*_{i}=*E*_{i}*ϵ*(*z*) for the longitudinal stress in the *i*th layer. Before rescaling, the tension is then given by [44]

Rescaling the constitutive relationships for the tension and the bending moment then yields

The Kirchhoff constraint *η*=0 (i.e. shear undeformability) results in the kinematic relations

When compared to the equations for a single-layer microbeam, coupling effects (captured by *k*_{12} multiplied by the square of the bending curvature) give rise to quadratic terms. These are due to the loss of elastic symmetry of the cross section for the bilayer geometry. This quadratic nonlinearity can mitigate the hardening effects of the geometric nonlinearity associated with midplane stretching [36].

Throughout this paper, we use parameter values *E*_{2}/*E*_{1}=0.312 and *ϱ*_{2}/*ϱ*_{1}=5.61 for the bilayer beam. Layers in this bilayer beam have identical widths, i.e. *A*_{2}/*A*_{1}=*h*_{2}/*h*_{1}.

## 3. Perturbation analysis for single-layer microbeams

To enable verification of the computational scheme used in §4 to analyse the forced response of a bilayer microbeam, we consider here the application of the multiple-scales perturbation method to the free response of a single-layer microbeam. In particular, we extend the analysis from Dankowicz & Lacarbonara [45] to second order in a small parameter, in order to capture higher order nonlinear effects. To this end, introduce a small non-dimensional parameter *ε*≪1, such that *α*^{2}↦*ε*^{−1}*α*^{2} (high beam slenderness), *v*↦*εv*, and *u*↦*ε*^{2}*u*. It follows that *ε*)

### (a) Method of multiple scales

We seek an oscillatory response of the two coupled PDEs. To this end, suppose
*t*_{i}=*ε*^{i}*t*, such that
*v*. Furthermore, let

We proceed to substitute these expansions for the unknown functions and the differential operators into the pair of coupled PDEs.

*Order* *ε*^{0}. To this order, we obtain
*n*th eigenfunction of the clamped–clamped Euler–Bernoulli beam corresponding to the natural frequency *ω*_{n}, i.e.
*n*th mode of the microbeam, we assume the generating solution of equation (3.7) as

Using the fixed–fixed boundary conditions

*Order* *ε*^{1}. To this order, we obtain
*q*=3, *q*=1 and *q*=1, respectively. We make the choices of *β*_{4} to 0.

We now have
*Order* *ε*^{2}. To this order, we obtain

### (b) Reconstitution

We express the complex-valued amplitude *A*_{n} in terms of real amplitude *a*_{n} and phase *φ*_{n} as *t*_{k}=*ε*^{k}*t*, it follows that the frequency of oscillation is given by
*α*=1, numerical evaluation of the integrals in equation (3.45) yields
*ε*, but positive at second order for the first mechanical mode of vibration. Furthermore, we see a hardening dependence on the response amplitude to first order in *ε*, but a softening contribution at second order.

## 4. Computational method

We proceed to develop a computational scheme for finding and continuing approximate solutions of equations (2.21) and (2.22) for clamped–clamped boundary conditions. To this end, proper discretizations of both spatial and temporal domains are required.

### (a) Mixed form

Before discretization, similar to our analysis of the scalar beam in [47], we proceed to rewrite equations (2.21) and (2.22) in the following mixed form:

### (b) Spatial and temporal discretization

The weak form of the governing equations is obtained by multiplying each of equations (4.1)–(4.5) by the corresponding test function *w*_{1}, *w*_{2}, *w*_{3}, *w*_{4} and *w*_{5}, and integrating over *Ξ*=[0,1]. For both trial and test functions, we choose continuously differentiable, piecewise-polynomial functions in *H*^{2} with arbitrary base point values on a spatial mesh given by the uniform partition 0<*s*_{1}<⋯<*s*_{M+1}=1 of the domain *Ξ*, and the uniform subpartition 0=*ξ*_{1}<⋯<*ξ*_{p+1}=1 of the *j*th interval *Ξ*_{j} in terms of the local variable *ξ*:=(*s*−*s*_{j})/(*s*_{j+1}−*s*_{j}). On *Ξ*_{j}, we express the unknown functions and test functions as
*w*_{1}, *w*_{2}, *w*_{5} vanish on the boundaries for the clamped–clamped beam. This choice results in a non-degenerate set of differential-algebraic equations of index-2. Note that the continuity conditions for *u*_{3}, *u*_{4} and *u*_{5} imply the compatibility conditions for the rotations of the cross sections, the bending moment and the shear force, respectively [48,49].

For the temporal discretization, we use a collocation method in terms of a continuous, piecewise polynomial function of degree *m* across *N* intervals, and impose that the residuals vanish at Radau points [47] along with the continuity conditions at the mesh points (cf. ch. 6 of [50]).

### (c) Closure conditions

Using parameter continuation, we investigate two types of dynamic behaviour of the beam, namely free (and undamped) and forced periodic vibrations. In each case, we append closure conditions on the discretized equations to ensure the local existence of a one-parameter family of periodic solutions near an initial solution guess. We use algorithms for covering implicitly defined manifolds in the computational continuation core [50] to successively impose projection conditions that generate a sequence of unique periodic responses along such a family.

The free vibration response gives the backbone of the forced frequency response. The simple imposition of periodicity results in a singular formulation, as there exist infinitely many periodic solutions for a given set of model parameter values, differing in amplitude. As discussed in [51], this issue can be resolved by introducing an additional unknown, *κ*, corresponding to artificial damping, so that the one-parameter family of periodic solutions is achieved locally only for *κ*=0.

In the case of forced vibrations, we perform continuation of a family of periodic responses to a harmonic excitation with *f*_{u}=0 and *Λ* is the amplitude of the excitation and *ω* is the corresponding frequency. To this end, similar to our analysis in [47], an autonomous pair of equations governing the time evolution of the excitation are appended to the discretized differential equations.

Finally, to resolve the degeneracy associated with arbitrary phase shifts along a periodic solution of an autonomous system, we treat the period *T* as unknown, rescale time by *T* and add the condition *t*=1.

### (d) Model verification

In this section, we consider the free undamped response of the single-layer beam and compare the results obtained from the numerical approach with those predicted by the second-order perturbation approximation in §3. We begin by considering the amplitude-independent corrections to the linear resonance frequency for the first mechanical mode, obtained in the limit of zero response amplitude. Table 1 shows sample results for three different values of *ε*, where *ε* becomes smaller, as predicted by the perturbation analysis. Indeed, the ratio *ε*.

Figure 2 shows the backbone curve for *ε*=10^{−3} using the perturbation method, along with that obtained using the numerical method. We observe close agreement with the predicted nonlinearity at moderate response amplitudes. The first-order perturbation (red line) predicts a higher bending of the backbone that arises from an overestimation of the hardening effects. This is mitigated by softening contributions captured at higher order, as shown by the green line.

## 5. Parametric study of the dynamic range

The geometric nonlinearities due to large-amplitude deflections result in three coexisting responses to excitations within certain frequency bandwidths. Although this nonlinear behaviour might be desirable in some applications, a linear response is often sought. Therefore, the larger the excitation amplitude at which there exists a single periodic solution in the frequency response, the larger the range of desired performance of the nano/microbeam resonator. The amplitude associated with the transition to nonlinear behaviour is called the *critical amplitude* and is denoted by *a*_{c}. In nanoresonators, the range of performance is often further limited by the fact that the noise-to-signal ratios are high. The difference between *a*_{c} and the noise level is called the *dynamic range* of the resonator [33,37,38,52]. Given a fixed noise level, both the critical amplitude and the dynamic range undergo identical variations. Throughout this paper, our focus is on optimizing the upper point of the dynamic range (i.e. the critical amplitude). The analysis is not concerned with changes to the dynamic range that result from design-induced variations in the noise level. Such a treatment is beyond the scope of this study.

### (a) An approximate cusp point

Using numerical continuation and the closure conditions described in §4, we obtain the frequency–response curve, a portion of which is shown in figure 3*a*, and for which the excitation amplitude is close to the onset of nonlinear behaviour. Here, *h*_{1}/ℓ=0.015, *h*_{2}=0, *Λ*=0.115, *c*_{u}=*c*_{v}=2*ζω*_{n}, *M*=12, *p*=3, *N*=10 and *m*=5, where *ω*_{n} is the natural frequency and *ζ*=0.02. This corresponds to the nonlinear regime, where more than one stable solution exists for a (relatively narrow) range of excitation frequencies, limited by two fold points, each of which is associated with a frequency and a particular periodic response. By reducing the excitation amplitude, this range of frequencies shrinks, until the two fold points merge at a cusp point, associated with a unique periodic response for some combination of excitation amplitude and frequency. In parameter space, this cusp point may be used as a threshold that distinguishes between linear and nonlinear regimes of operation.

To approximately locate the cusp point, we proceed to define four instances of the boundary-value problem with initial solution guesses corresponding to the red points in figure 3*a*. Here, *v*_{l} and *v*_{r} denote the amplitude differences for the left and right pair, respectively. We obtain a composite continuation problem by imposing additional constraints on the four instances of the forced response boundary-value problem. First, we require that each pair of solutions have identical frequencies throughout the continuation. Then, we fix *Λ* so that the underlying frequency–response curve is not changing. We perform the continuation by varying *ω*_{r}, while *ω*_{l} is kept fixed, until *v*_{r}=Δ_{v}, where Δ_{v} is a relatively small number indicating the prescribed tolerance. The same procedure is repeated for *ω*_{l}, while *ω*_{r} is fixed, until *v*_{l}=*v*_{r}. Next, for fixed *v*_{l} and *v*_{r}, we define *ω*_{δ}=*ω*_{r}−*ω*_{l} and perform continuation until *ω*_{δ}=0 (figure 3*b*), allowing the amplitude of the excitation to vary. The four periodic orbits and the corresponding excitation frequency and the excitation amplitude obtained at the end of continuation approximate the cusp point. Notably, if Δ_{v} approaches zero, the four periodic orbits merge at the cusp point. In this limit, however, the problem becomes ill-conditioned and, ultimately, singular. Here, to avoid this singularity, by accepting a slight overestimation of the dynamic range, we use Δ_{v}=10^{−4} and let the amplitude corresponding to the upper point represent *a*_{c}.

To allow for continuation of the approximate cusp point, we substitute a composite continuation problem with three instances of the forced response boundary-value problem with initial solution guesses given by the three distinct periodic solutions found above. Here, all three responses are required to have identical frequency and *v*_{l}=*v*_{r}=Δ_{v} throughout continuation.

### (b) Numerical study

We proceed to study the effects of the beam geometry on the critical amplitude *a*_{c}. Figure 4 shows that *a*_{c} increases approximately linearly with respect to *h*_{1}/ℓ for a single-layer microbeam. This increase is likely associated with an increased softening contribution from rotary inertia due to the change in the aspect ratio of the cross section.

We choose parameter value *h*_{1}/ℓ=0.01 and add a uniform layer on top of the beam’s substrate. In addition, we consider two cases in which the added layer is non-uniform and has step changes in its thickness at a distance ℓ/3 from each edge. The case of step changes corresponds, for example, to actual designs where piezo layers are added for actuation or sensing purposes. In the first case, the added layer covers the outside sections (i.e. the middle section remains a single layer) while in the second case, only the middle section is bilayered. We refer to these two cases by *YNY* and *NYN*, respectively, and the case where the top layer covers the entire span by *YYY*. Figure 5 shows the corresponding relationship between *h*_{2}/*h*_{1} and *a*_{c}, for *h*_{1}/ℓ=0.01. Within the range shown in the figure, the optimal parameter choice, corresponding to maximal *a*_{c}, occurs when *h*_{2}=0.5*h*_{1} in the *YNY* case. Indeed, the value of *a*_{c} corresponding to this configuration is consistently higher than for either of the other two configurations. Interestingly, the uniform beam has smaller dynamic range than either of the other two configurations throughout this parameter range.

To explore this behaviour further, we show in figure 6 the graphs of *u*_{2}, *u*_{3}, *u*_{4} and *u*_{5} for *h*_{1}/ℓ=0.01 and *h*_{2}/*h*_{1}=0.5, at the moment in time when the midpoint experiences maximum lateral deflection. In each panel, the black curve corresponds to *h*_{2}=0, i.e. the single-layer beam. Notably, the behaviour of the *NYN* case is close to that of the single-layer beam in the outside sections and transitions to that of the uniform beam in the middle section. Similarly, the behaviour of the *YNY* case is close to that of *YYY* case in the outside section and transitions to the single-layer beam behaviour in the middle section. It appears that it is only in the case where a layer is added to the two ends of the beam, that the beam *softens* and exhibits higher deflections compared with the single-layer beam.

Oscillations are evident in the graph of *u*_{5}, and to a lesser extent *u*_{4}, in the cases of discontinuous changes in the beam axial and bending stiffnesses. To accurately describe the associated discontinuous changes in elongation and curvature requires mesh refinement in the vicinity of the discontinuities. Such local effects, while propagated along the entire domain in *u*_{5}, tend to become less significant when translated to lower orders, i.e. *u*_{2}, *u*_{3} and *u*_{4}. This is illustrated in figure 7, where the graphs of *u*_{2}, *u*_{3}, *u*_{4} and *u*_{5} are shown in the *YNY* case for different number of equal-sized finite elements. We observe that although the oscillatory behaviour of *u*_{5} is directly related to the number of elements, convergence to a nonoscillatory behaviour is seen as the number of elements increases. Further evidence of convergence is shown in table 2, where the values of the non-dimensional stored energy equation (2.18) are obtained for each beam configuration and each discretization order.

Finally, we proceed to investigate the dependence of the dynamic range on the length of the top layer. To this end, we again consider two cases: *NYN* and *YNY*. However, the span of the beam is no longer divided into three equal segments. Let *d* denote the accumulated length of the top layer. For the *NYN* case (figure 8*a*), the dynamic range has a global maximum when *d*=0. For *h*_{2}/*h*_{1}=0.5, there are two local maxima; one at *d*≈0.4 and the other one at *d*=1, corresponding to a uniform bilayer beam. In the case when *h*_{2}/*h*_{1}=0.2, there is a global minimum at *d*≈0.7.

Figure 8*b* shows the results for the *YNY* case. Here, a global maximum is found away from *d*=0 and is significantly higher when *h*_{2}/*h*_{1}=0.5 than *h*_{2}/*h*_{1}=0.2. We observe a 10–25% variation in the dynamic range depending on the specific placement of the top layer, and further note that the *YNY* case always gives a larger dynamic range than the *NYN* case.

## 6. Conclusion

The scheme developed in this paper has a broad range of applications in the design of N/MEMS devices. In fact, beam structures are present in many of such systems due to their relatively straightforward fabrication process, robust performance, as well as their versatility as a functional part at small scales. In all of these applications, high accuracy predictions along with precise measurements of system’s behaviour are vital in achieving a reliable device. Consequently, performance in linear regimes has been an essential element of their design. These applications encompass almost all N/MEMS devices including actuators [53], oscillators [54], force sensors [55], filters [56], magnetic resonance force microscopy [57], ultrasensitive mass detection [58] and single spin detection [59]. Despite such ubiquity, there is no general computational scheme that can be employed in nonlinear tuning of N/MEMS devices.

As suggested by the analysis, care should be taken to account for the coupling between longitudinal and lateral dynamics, especially in instances where large deflection leads to significant geometrically nonlinear effects. As further argued and demonstrated in this paper, computational tools should be verified against analytical results, as available, prior to use in design optimization.

The application of the methodology to a bilayered microbeam is motivated by actual device designs that rely on a second layer, integrated on top of the resonator substrate, to act as an actuating or sensing element. As shown in the numerical results, the geometry of the added layer may affect the dynamic range of the resonator by a significant amount. The formulation of an inverse problem associated with the near-cusp behaviour of the frequency–response curve provides a rigorous foundation for analysing the parameter sensitivity of the dynamic range using numerical continuation.

## Author' contributions

W.L. and H.D. conceived the study, derived the microbeam model and designed the perturbation analysis. M.S. implemented the perturbation and continuation analysis. M.S. further conceived the dynamic range study. All authors participated in writing the manuscript and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This material is based on work supported by the National Science Foundation under grant nos. 1016467 and 0855787. The work of Walter Lacarbonara was partially supported by the Italian Ministry of Education, University and Scientific Research (2010-2011 PRIN grant no. 2010BFXRHS-002).

- Received December 16, 2014.
- Accepted June 1, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.