Nonlinear phenomena such as internal resonances have significant potential applications in micro electro mechanical systems (MEMS) for increasing the sensitivity of biological and chemical sensors and signal processing elements in circuits. While several theoretical systems are known which exhibit 1:2 or 1:3 internal resonances, designing systems that have the desired properties required for internal resonance and that are physically realizable as MEMS devices is a significant challenge. Traditionally, the design process for obtaining resonant structures exhibiting an internal resonance has relied heavily on the designer's prior knowledge and experience. However, with advances in computing power and topology optimization techniques, it should be possible to synthesize structures with the required nonlinear properties (such as having modal interactions) computationally. In this work, a preliminary work on computer based synthesis of structures consisting of beams for desired internal resonance is presented. The linear structural design is accomplished by a Finite Element Method (FEM) formulation implemented in Matlab to start with a base structure and iteratively modify it to obtain a structure with the desired properties. Possible design criteria are having the first two natural frequencies of the structure in some required ratio (such as 1:2 or 1:3). Once a topology of the structure is achieved that meets the desired criterion, the linear mode shapes of the structure can be extracted from the finite element analysis, and a more complete Lagrangian formulation of the nonlinear elastic structure can be used to develop a nonlinear two-mode model of the structure. The reduced-order model is expected to capture the appropriate resonant dynamics associated with modal interactions between the two modes. The nonlinear response of the structure can be obtained by application of perturbation methods such as averaging on the two-mode model. Many candidate structures are synthesized that meet the desired modal frequency criterion and their nonlinear responses are compared.

## Introduction

In the past decade, there has been an increased emphasis on using nonlinear dynamics in the design of MEMS devices [1–3]. In particular, the phenomenon of internal resonance has recently attracted attention due to the several advantages it can provide for MEMS based sensing [4,5]. These advantages include extreme sensitivity to perturbations in the measured quantity such as mass and ease of measurement [5]. Internal resonance is a nonlinear phenomenon in which transfer of energy can take place between modes depending on the type of nonlinearity present; for example if a system has quadratic nonlinearities and has two of its modal frequencies in ratio of 1:2, then it can exhibit 1:2 internal resonance if the excitation amplitude is above some threshold [6–8]. However, the design of structures for internal resonance is somewhat constrained by the designer's experience. Topology optimization is a popular synthesis technique already in use for design of compliant mechanisms [9]. It has also been used in MEMS design for example to get devices with particular frequency characteristics [10,11]. One of the biggest hurdles in design of structures for internal resonance is to get their natural frequencies in some ratio (1:2 or 1:3) and to then see if the appropriate modes are nonlinearly coupled to exhibit nonlinear interactions.

This work describes a preliminary computational synthesis effort based on FEM to achieve natural frequencies of a structure in some desired ratio that helps in the initial design stages of a MEMS device which uses internal resonance as it is working principle. While in theory, energy transfer can occur between any two different modes of a nonlinear structure, the present work is restricted to considering structures with energy transfer between the first two modes. Generally, the higher a mode's natural frequency, the larger the energy required to excite that particular mode in steady state. Therefore, the lowest two modes of a structure are more easily excited.

The structures obtained by the process outlined in this work are all composed of flexible beams that are rigidly connected orthogonally and undergo planar vibrations. A hierarchical optimization procedure is adopted to obtain the structures having their lowest two natural frequencies in some desired ratio. Once a candidate structure has been identified, its mode shapes can also be extracted as they are available through the FEM. The section on Linear Structure Synthesis describes the optimization procedure and presents some examples of final structures having the desired properties. The section on Nonlinear Frequency Response develops the response of the systems synthesized through the optimization process. A two-mode model of the structure is constructed by expressing the transverse deflection of the beams as a linear combination of the lowest two modes. Only a two-mode model is considered because it is assumed that energy transfer occurs only between these two modes. More specifically, the second mode is directly excited using external excitation and it is hoped that it in turn excites the first mode due to nonlinearities present in the structure. While in general the nonlinear response of the structure may consist of several other modes, it is expected that in the presence of damping, modes with neither a direct excitation nor an internal resonance will have their amplitudes diminish over time [12]. Therefore, a two-mode model can be assumed to provide a fairly accurate representation of the system's nonlinear response. This two-mode approximation then leads to a two degree-of-freedom reduced-order model that can be used to obtain equations for the slow time amplitudes [6,12]. These slow-time amplitude equations can be solved for obtaining the final nonlinear frequency response of the structure. Finally, the Summary and Conclusions section discusses potential applications and improvements in the methodology.

## Linear Structure Synthesis

As mentioned in the introduction, the structures which are initially synthesized by this method consist of thin beams undergoing planar vibrations. The natural frequencies and mode shapes of these structures are obtained using the FEM analysis with planar two-node beam elements [13]. Each node has three degrees of freedom namely longitudinal, transverse, and rotational. The transverse and rotational degrees of freedom are interpolated using cubic polynomials and the longitudinal degree of freedom is interpolated using linear polynomials.

The major purpose of the linear structure synthesis is to obtain a structure which has its lowest two natural frequencies in some desired ratio (such as 1:2 or 1:3). The inputs to the optimization process are:

- (1)
Desired natural frequency ratio (1:2 or 1:3).

- (2)
Desired frequency range (a lower and an upper bound on the frequencies).

where Γ_{1} and Γ_{2} are the lower and upper bounds of the desired frequency range. Based on these, the optimization problem can be formulated as

A hierarchical optimization process is used to synthesize the structures. It is assumed that the entire structure is built out of the same material and that all beams have the same material properties such as density and elastic modulus. For results in this work, it was assumed that the beams were made of Silicon (Si) with Young's Modulus 185 GPa and density equal to 2329 kg/m^{3} [14]. In the first-level optimization, a base structure is chosen. The base structure may be a simple cantilever or may have several beams connected together in some configuration. Four of the structures which were used as base structures in the optimization procedure are shown in Fig. 1. All the beams in the base structure have the same length and the same cross-section. The base structure is then modified by adding a single beam at a time which has the same length, cross-section, and material properties as the beams of the base structure. In the present work, beams can only be added orthogonally to existing beams and only at their end points. Normally, there is more than one option available to add a single beam to a base structure as a beam can be added at multiple points and in multiple configurations. Adding a beam in each one of these possible options leads to a new structure. The objective function value *c*(ω) is calculated for each of these new structures and the structure with the least objective function value is chosen for further improvement, while assuring that this least objective function value is less than the objective function value of the original structure. If the *c*(ω) value for the new structure is less than the *c*(ω) value for the original structure, then again all possible options for adding a single beam to the new structure are evaluated and the best option is chosen. This process is continued until no further reduction in the objective function value is possible by adding a beam of the same length and cross-section. The structure is now subjected to the second-level optimization.

In the second-level optimization, the lengths of individual beams of the structure which was finally obtained from the first-level optimization are optimized using a gradient based optimization method (more specifically the method of steepest descent [15]) to obtain further reduction in the objective function (*c*(ω)) value. All beam lengths are optimized simultaneously. In most cases, during the course of this work, as will be shown in an example presented later, structures with the same number of beams connected in the same manner but with unequal beam lengths resulted in lower objective function (*c*(ω)) values than structures with equal beam lengths. The second-level optimization is stopped when either the objective function gradient becomes sufficiently small or the objective function value falls below the specified tolerance limit. It is thus possible to view the first-level optimization as the shape optimization of the structure and the second-level optimization as the size optimization of the structure. After the beam lengths have been optimized, the structure is subjected to the third and final-level optimization.

In the third-level optimization, the structure is made to satisfy the bounds constraint (Eq. (2)). This is achieved by changing the dimensions of the entire structure, such as the beam lengths and cross-sectional properties, by the same amount, i.e., by scaling the dimensions of the entire structure appropriately so that the second natural frequency lies within the specified limits. Note that scaling the dimensions such as beam lengths and cross-sectional properties of the entire structure by the same amount changes individual natural frequencies but does not change their ratio.

Overall, this hierarchical optimization procedure can be summarized as follows:

- (1)
The natural frequencies of the base structure are calculated using FEM. In this step, the entire structure is divided into a number of beam elements while keeping in mind that the length of each element be at least ten times larger than its width or depth so as to minimize deviation from Euler–Bernoulli beam behavior. Let

*M*and*K*be the final assembled mass and stiffness matrices of the base structure, respectively; then the eigenvalues of the matrix*M*^{−1}*K*are the squares of the natural frequencies of the structure. Knowing the natural frequencies the objective function*c*(ω) can be evaluated. - (2)
The base structure is modified by adding another beam to it. This additional beam is of the same length and cross-section as the beams of the base structure and can be added only at the unclamped start or end nodes of individual beams. Also, this beam can be added only in orientations of ±90 deg or 180 deg with respect to existing beams in the structure. Often there are several such nodes and configurations possible for the additional beam and for each of these possibilities the objective function

*c*(ω) is evaluated to determine the structure to be used for the next iteration. - (3)
If the incorporation of any additional beam in any possible configuration does not lead to an objective function value less than the one for the present structure, the present structure is taken as the optimum shape. The lengths of the individual beams are now optimized by the gradient based method of steepest descent using a program which the authors have developed. The optimization is stopped when the numerical gradient gets sufficiently small or the objective function value is below a specified tolerance limit.

- (4)
Finally, if the second natural frequency ω

_{2}does not lie within the specified limits, the dimensions of the entire structure (lengths and cross-sections) can be scaled by the same value to shift the natural frequencies within the required limits.

A flow chart depicting the important steps in the optimization procedure is shown in Fig. 2.

### An Example.

For illustration purposes, let it be assumed that a structure needs to be designed for 1:2 resonance starting with a cantilever beam as the base structure. Furthermore, let the desired second natural frequency (ω_{2}) be between 3 × 10^{5} and 3.5 × 10^{5 }rad/s (between about 50–55 kHz). Thus, as can be seen in Fig. 3(a), the process starts with a simple cantilever beam. The subsequent steps of the first-level optimization are shown in Figs. 3(b)–3(f). In Fig. 3(b), the three possible configurations for the additional beam are shown with dotted lines. Figure 3(c) shows the optimal structure with two beams (solid lines), and the possible configurations of an additional beam to be added. Continuing in this manner, the structure in Fig. 3(e) is obtained. After the structure in Fig. 3(e) is obtained, none of the news additions of beams shown with dotted lines leads to a further reduction in the objective function. Thus, the final structure obtained from the first level optimization is given by Fig. 3(f). After this, the lengths of the beams are optimized using the method of steepest descent to further improve the objective function in the second-level optimization.

The objective function value (Eq. (3)) for the structure in Fig. 4(a) is 0.0063 while the corresponding value for the structure in Fig. 4(b) is 8.003 × 10^{−11}. As can be observed from Table 1, slight changes to individual beam lengths leads to a significant improvement in the objective function. The final results obtained from the optimization process are presented in Table 2.

### Verification and Mode Shapes.

In order to verify the FEM program used for structural optimization implemented in this preliminary work, the natural frequencies of the structure shown in Fig. 5 were computed using the present formulation as well as the commercially available FEM software ANSYS. This is the same structure as shown in Fig. 4(a). In both cases, the structure was divided into the same type of elements (2-node beam element) and the same number of nodes and elements were used. As can be seen from Table 3, the percentage errors in natural frequency are negligible.

*M*

^{−1}

*K*mentioned earlier. For ease of calculations, the deflections of individual beams are represented as polynomials; thus, for the

*i*th beam of a structure let its arc length be denoted by

*s*and its transverse deflection in the

_{i}*j*th mode shape be denoted by

*v*. The relationship between the transverse deflection and the arc length is assumed to be of the form

_{ij}where *N* is the degree of the polynomial. The coefficients *b _{ijk}* are obtained by curve fitting the deflections obtained using FEM. This curve fitting is done using the “polyfit” function of Matlab, which fits a specified degree polynomial to a set of data in the least-square sense.

For verification of the mode shapes, the lowest two natural frequencies and the corresponding mode shapes of the structure shown in Fig. 4(b) were obtained using Hamilton's Principle [5] and compared with the results obtained from FEM and curve fitting. For a more exact description of the application of Hamilton's Principle to the structure shown in Fig. 4(b), please refer to the section on the Nonlinear Frequency Response. Table 4 lists the lowest two natural frequencies obtained by the two methods. The corresponding modes of the structure in Fig. 4(b) are shown in Fig. 6.

The verification of mode shapes presents some difficulties as there as several ways to compare the mode shapes. In this work, the tip deflection of one of the beams was made equal in mode shapes obtained from both methods (FEM and Hamilton's Principle). Figure 7 presents the comparison between the mode shapes obtained using FEM and curve fitting and Hamilton's Principle; as can be seen, the match is excellent and therefore the FEM program is assumed to give fairly accurate natural frequencies and mode shapes. To get a measure of the relative error between the two mode shapes obtained using the two methods, the following method was adopted; let *v _{ij}* be the transverse deflection of the

*i*th beam for the

*j*th mode obtained using FEM and curve fitting and let

*w*be the same quantity obtained using Hamilton's Principle. Then a modal error quantity

_{ij}*E*may be defined by

_{j}where *BN* is the number of beams in the structure. For the structure shown in Fig. 4(b), the error for the lowest mode (*E*_{1}) was 1.1407 × 10^{−22} m and the error for the second mode (*E*_{2}) was 2.8054 × 10^{−22} m, with the beam lengths in the structure being of the order of 10^{−4} m.

### Results for Linear Synthesis.

The final structures for 1:2 internal resonance obtained by starting from the base structures in Fig. 1 are given in Fig. 8. For later reference, let the structures in Fig. 8 be called as structure 1, 2, 3, and 4, respectively. Table 5 presents the important frequency properties of these structures.

The same base structures can also be used for design of structures with other requirements, e.g., for 1:3 internal resonance. The structures 5, 6, 7, and 8 shown in Fig. 9 have been obtained by starting from the base structures shown in Fig. 1. Table 6 presents the important frequency properties of the structures 5, 6, 7, and 8.

## Nonlinear Frequency Response Prediction

Consider the structure in Fig. 8(a); the nonlinear responses of the other structures shown in Fig. 8 can be obtained in a similar manner. Figure 10 gives a more complete schematic of structure 1 (see Table 7).

*T*and the strain energy

*U*of the structure can be written as [8,16].

*v*is the transverse deflection,

_{i}*u*is the axial deflection, and ρ

_{i}*is the linear mass density of the*

_{i}*i*th beam, and

where *E* is the material's Youngs Modulus, *I _{i}* are the area moments of inertia, and

*A*are the cross sectional areas of the beams.

_{i}*A*

_{1}and

*A*

_{2}are the modal amplitudes that are assumed to be only functions of time, and ϕ

_{i}_{1}and ϕ

_{i}_{2}are the first and second mode shapes of the

*i*th beam, respectively. As already discussed, the mode shapes can be obtained as polynomials. One particular point which must be mentioned here is that the mode shape polynomials are obtained by curve fitting on the scaled mode shapes. The scaling is assumed to be 10% of the first beam length; therefore, the mode shape polynomials have a dimension associated with them and consequently the modal amplitudes,

*A*

_{1}and

*A*

_{2}are nondimensional. Using the inextensibility constraint (Eq. (8)), the axial displacement can be obtained as [7]

where *p _{i}* and

*q*are amplitude components which vary on a slow time scale τ = ε

_{i}*t*, as defined in Refs. [6,7].

*L*) of the structure can now be obtained by substituting Eqs. (9), (10), (12), and (15) and their derivatives into the expressions of the kinetic and strain energies of the structure defined in Eqs. (6) and (7). This Lagrangian is then averaged over the period of the oscillation $4\pi \Omega $. The amplitudes

*p*and

_{i}s*q*and their derivatives w.r.t. τ, are treated as constants in this averaging as they are functions of the slow time scale. In addition, only terms up to

_{i}s*O*(ε

^{3}) are retained as 1:2 internal resonance is primarily a result of coupling due to quadratic nonlinearities. The averaged Lagrangian of the system is defined by

*p*and

_{i}s*q*:

_{i}swhere a (′) now denotes a derivative with respect to the slow time τ.

Note that the equations for the slow-time amplitudes have been modified to include modal damping terms (ζ_{1} and ζ_{2}). The modal damping coefficients (ζ_{1} and ζ_{2}) were set equal to 0.1 arbitrarily as no experiments were conducted to estimate the real damping coefficients. Equations (17*a*)–(17*d*) are similar to the equations for modal amplitudes obtained in Refs. [6,7,12] for the case of 1:2 internal resonance in a system with two degrees of freedom. These equations can be solved for steady-state solutions to give single-mode (only second modal amplitude is nonzero) and coupled-mode solutions (both first and second mode amplitudes are nonzero). Figure 11 illustrates some important aspects of the nonlinear dynamic response of structure 1. In a linear system, when a particular mode is excited at resonance, it acquires a large amplitude which theoretically, in the absence of damping, can even go up to infinity; correspondingly, the contributions of the other modes to the system response are zero or minimal at that particular modal frequency. This is what is observed in the second mode's linear response in Fig. 11, where it can be clearly seen that the second mode has a large amplitude when excited at its natural frequency. However, in the presence of nonlinearities and when the frequencies of the first two modes of the system are in a specific ratio, an energy transfer takes place between the second and first modes which again can be observed from Fig. 11, wherein the second mode's amplitude in the nonlinear response is attenuated (as compared to the linear response) and the first mode appears with a finite, nontrivial amplitude even when the system is being excited near its second natural frequency. The coupled-mode solution is of interest here as it represents the true solution with modal coupling resulting from internal resonance. Despite the structure being excited close to its second natural frequency, a nonzero response of the first mode at almost half the excitation frequency is obtained due to quadratic nonlinearities in the system. Figure 12 shows the magnitude of the first-mode amplitude (defined by $a1=p12+q12$) obtained for structures 1,2,3, and 4 shown in Fig. 8. The responses are shown all together in Fig. 13 for a quantitative comparison as well. For another quantitative comparison of the candidate structures developed, the point of maximum deflection for the first mode for each of the structures 1, 2, 3, and 4 is taken and its total deflection is computed using Eq. (9). Figure 14 shows the location of the point of maximum deflection for the first mode for structures 1, 2, 3, and 4. As can be observed from Fig. 14, the point of maximum deflection for the first mode occurs at the end of one of the constituent beams and thus its deflection can also be referred to as the tip deflection. These tip deflections will vary periodically with time with a constant amplitude during steady state. Figure 15 shows the comparison between the tip amplitudes for the structures 1, 2, 3, and 4 when they are subjected to the same level of base excitation. It can be observed from Fig. 15 that different structures have different magnitudes of deflection for the same level of excitation and thus from an application standpoint, some structures may be more appropriate than others.

## Discussion and Future Work

A computational method based on a hierarchical optimization procedure has been described to aid the design process for MEMS devices which utilize nonlinear phenomena such as internal resonance as their operating principle. Several examples of structures are presented which are possible candidates for exhibiting 1:2 or 1:3 internal resonance in two of their linear structural modes. Although not presented, the method can also be used in the design of structures which might potentially exhibit 1:4 or other higher-order internal resonances. For the case of 1:2 internal resonance, the nonlinear frequency response of several structures developed using the classical approach of the two-mode models and averaged Lagrangian is presented. As is evident, using the method proposed, a large number of structures displaying the same kind of qualitative nonlinear behavior can be developed. This then provides a large database from which to pick the desired design based on their quantitative considerations; for example, as can be observed from Fig. 12 and more clearly from Fig. 15, for the same level of excitation, structure 4 provides much larger tip amplitude than other structures thus showing most promise among the four structures for applications where high response amplitudes are required. Note that this may be partly due to relatively larger dimensions of structure 4. Furthermore, note that in the case of other types of internal resonances, the nonlinear Lagrangian may need to take into account extensible beam models such as stretching of the centerline which can be certainly incorporated without much difficulty. It should also be pointed out that including more than just the two internally resonant modes for the reduced order model can lend further to the analysis, particularly in capturing the softening effects for the nonlinear response of the structures. However, this will not change the overarching conclusions or comparisons of responses for the structures presented in this work.

The current work can be enhanced in different directions. It could be extended to synthesize nonplanar structures. In the optimization procedure itself, at the first level, provision could be made to add beams in any configuration to the base structure. Moreover, some flexibility could be built-in with regards to the dimensions of the beam being added. In the second-level optimization, techniques which optimize all beam lengths as well as beam cross-sections concurrently rather than just the beam lengths could be implemented; however, beams of the same cross-section give advantages in terms of general ease of fabrication of the final structure. Additionally, elements other than beams such as plates can be incorporated to further enhance the number of options available for a particular kind of device. Nonstructural elements such as electrostatic actuator drives common in many MEMS devices can also be made part of device models to perform better device simulations and optimization. Finally, multilevel optimization can be pursued to also include the amplitude of nonlinear response as an objective. Some of these avenues are currently being pursued and the results will be presented in forthcoming articles.

## Acknowledgment

The authors would like to acknowledge the financial support provided by Fredrick N. Andrews Fellowship and the Alpha P. Jamison Professorship funds at Purdue University. The authors also thank the reviewers for a careful review and constructive suggestions.