Computational Synthesis for Nonlinear Dynamics Based Design of Planar Resonant Structures

Tripathi, Astitva; Bajaj, Anil K.

doi:10.1115/1.4024845

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).

Analytically, the above requirements can be specified by the relations

\frac{ω_{2}}{ω_{1}} = η

(1)

where ω₁ and ω₂ are the first and second natural frequency of the structure, respectively, and η is the ratio specified in the input (η = 2 if the required ratio is 1:2), and

Γ_{1} \leq ω_{2} \leq Γ_{2}

(2)

where Γ₁ and Γ₂ are the lower and upper bounds of the desired frequency range. Based on these, the optimization problem can be formulated as

minimize

\begin{matrix} c (ω) = {(η - \frac{ω_{2}}{ω_{1}})}^{2} subject to Γ_{1} \leq ω_{2} \leq Γ_{2} \end{matrix}

(3)

where c(ω) is the objective function that has to be minimized.

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³ [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.

Fig. 1

View large Download slide

Some of the starting base structures for the optimization process (all dimensions are in meters)

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⁻¹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 ω₂ 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.

Fig. 2

View large Download slide

A flowchart of the optimization process

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 (ω₂) be between 3 × 10⁵ and 3.5 × 10⁵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.

Fig. 3

View large Download slide

Iterations of the optimization process to achieve 1:2 internal resonance. Actual structure is denoted by solid lines. Dotted lines denote possible configurations in which the extra beam may be attached (all dimensions are in meters).

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⁻¹¹. 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.

Fig. 4

View large Download slide

Local optimization of individual beam lengths to improve the structure (all dimensions are in meters)

Table 1

Beam lengths for the structure in Fig. 4, before and after application of the second-level optimization

Beam length	Before (Fig. 4(a))	After (Fig. 4(b))
L₁ (m)	1.0000 × 10⁻⁴	1.0158 × 10⁻⁴
L₂ (m)	1.0000 × 10⁻⁴	9.5351 × 10⁻⁵
L₃ (m)	1.0000 × 10⁻⁴	1.0403 × 10⁻⁴
L₄ (m)	1.0000 × 10⁻⁴	9.8892 × 10⁻⁵

Beam length	Before (Fig. 4(a))	After (Fig. 4(b))
L₁ (m)	1.0000 × 10⁻⁴	1.0158 × 10⁻⁴
L₂ (m)	1.0000 × 10⁻⁴	9.5351 × 10⁻⁵
L₃ (m)	1.0000 × 10⁻⁴	1.0403 × 10⁻⁴
L₄ (m)	1.0000 × 10⁻⁴	9.8892 × 10⁻⁵

Table 2

Natural frequencies for the structure in Fig. 4(b), obtained by the optimization process for 1:2 internal resonance with cantilever beam as the starting base structure

Variable	Value
ω₁ (rad/s)	1.6000 × 10⁵
ω₂ (rad/s)	3.2000 × 10⁵
$\frac{ω_{2}}{ω_{1}}$	2.0000

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.

Fig. 5

View large Download slide

Structure used for verification of the fem formulation with Ansys (all dimensions are in meters)

Table 3

Comparison of natural frequencies obtained by Ansys and Matlab based fem for the structure in Fig. 5

Data	Present program	Ansys
ω₁ (rad/s)	1.5737 × 10⁵	1.5738 × 10⁵
ω₂ (rad/s)	3.2723 × 10⁵	3.2769 × 10⁵
ω₃ (rad/s)	7.6552 × 10⁵	7.6717 × 10⁵
ω₄ (rad/s)	8.9701 × 10⁵	8.9981 × 10⁵

Data	Present program	Ansys
ω₁ (rad/s)	1.5737 × 10⁵	1.5738 × 10⁵
ω₂ (rad/s)	3.2723 × 10⁵	3.2769 × 10⁵
ω₃ (rad/s)	7.6552 × 10⁵	7.6717 × 10⁵
ω₄ (rad/s)	8.9701 × 10⁵	8.9981 × 10⁵

The formulation can also be used to obtain the mode shapes of the structure which are essential for constructing the two-mode reduced-order model of the nonlinear system. Figure 4(b) shows the structure which was obtained by the optimization procedure outlined in this paper for 1:2 internal resonance by starting with a cantilever as the base structure. The mode shapes are the eigenvectors of the matrix M⁻¹K mentioned earlier. For ease of calculations, the deflections of individual beams are represented as polynomials; thus, for the ith beam of a structure let its arc length be denoted by s_i and its transverse deflection in the jth mode shape be denoted by v_ij. The relationship between the transverse deflection and the arc length is assumed to be of the form

v_{ij} (s_{i}) = \sum_{k = 0}^{N} b_{ijk} s_{i}^{k}

(4)

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.

Fig. 6

View large Download slide

Lowest two mode shapes of the structure in Fig. 4(b) (all dimensions are in meters)

Table 4

Comparison of natural frequencies of the structure in Fig. 4(b) obtained by FEM and Hamilton's Principle

Variable	FEM	Hamilton's Principle
ω₁ (rad/s)	1.6000 × 10⁵	1.6000 × 10⁵
ω₂ (rad/s)	3.2000 × 10⁵	3.2000 × 10⁵

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 ith beam for the jth mode obtained using FEM and curve fitting and let w_ij be the same quantity obtained using Hamilton's Principle. Then a modal error quantity E_j may be defined by

Fig. 7

View large Download slide

Verification of the lowest two mode shapes obtained via FEM. There are two sets of lines representing mode shapes obtained using Hamilton's Principle and mode shapes obtained using fem (all dimensions are in meters).

E_{j} = \sqrt{\sum_{i = 1}^{BN} \int_{0}^{l_{i}} {(v_{ij} - w_{ij})}^{2} {ds}_{i}}

(5)

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₁) was 1.1407 × 10⁻²² m and the error for the second mode (E₂) was 2.8054 × 10⁻²² m, with the beam lengths in the structure being of the order of 10⁻⁴ 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.

Fig. 8

View large Download slide

Optimized structures obtained for 1:2 internal resonance by starting from the base structures shown in Fig. 1 (all dimensions are in meters)

Table 5

Frequency properties of structures with 1:2 internal resonance, as shown in Fig. 8

Structure	ω₁ (rad/s)	ω₂ (rad/s)	$ω_{2} / ω_{1}$
1	1.6000 × 10⁵	3.2000 × 10⁵	2.0000
2	1.6000 × 10⁵	3.1999 × 10⁵	2.0000
3	1.6000 × 10⁵	3.1999 × 10⁵	2.0000
4	1.6000 × 10⁵	3.1999 × 10⁵	2.0000

Structure	ω₁ (rad/s)	ω₂ (rad/s)	$ω_{2} / ω_{1}$
1	1.6000 × 10⁵	3.2000 × 10⁵	2.0000
2	1.6000 × 10⁵	3.1999 × 10⁵	2.0000
3	1.6000 × 10⁵	3.1999 × 10⁵	2.0000
4	1.6000 × 10⁵	3.1999 × 10⁵	2.0000

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.

Fig. 9

View large Download slide

Structures obtained for 1:3 internal resonance by starting from the base structures shown in Fig. 1 (all dimensions are in meters)

Table 6

Frequency properties of structures obtained for 1:3 internal resonance, as shown in Fig. 9

Structure	ω₁ (rad/s)	ω₂ (rad/s)	$ω_{2} / ω_{1}$
5	1.0666 × 10⁵	3.1999 × 10⁵	3.0001
6	1.0666 × 10⁵	3.1999 × 10⁵	3.0001
7	1.0666 × 10⁵	3.1999 × 10⁵	3.0001
8	1.0666 × 10⁵	3.1999 × 10⁵	3.0001

Structure	ω₁ (rad/s)	ω₂ (rad/s)	$ω_{2} / ω_{1}$
5	1.0666 × 10⁵	3.1999 × 10⁵	3.0001
6	1.0666 × 10⁵	3.1999 × 10⁵	3.0001
7	1.0666 × 10⁵	3.1999 × 10⁵	3.0001
8	1.0666 × 10⁵	3.1999 × 10⁵	3.0001

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).

Fig. 10

View large Download slide

A schematic diagram for the optimal structure 1 in Fig. 8 (all dimensions are in meters)

Table 7

Properties of the structure 1 shown in Fig. 10

Variable	Value
l₁ (m)	1.0158 × 10⁻⁴
l₂ (m)	9.5351 × 10⁻⁵
l₃ (m)	1.0403 × 10⁻⁴
l₄ (m)	9.8892 × 10⁻⁵
w (m)	1 × 10⁻⁶
h (m)	1 × 10⁻⁶
E (N/m²)	1.85 × 10¹¹
ρ (kg/m³)	2329

Structure 1 consists of four beams, each having the same cross-section and material properties. The structure is also subjected to a base excitation of

\overset{··}{w}

⁠. The kinetic energy T and the strain energy U of the structure can be written as [8,16].

\begin{matrix} T = \frac{1}{2} \int_{0}^{l_{1}} ρ_{1} [{({\overset{\cdot}{v}}_{1} + \overset{\cdot}{w})}^{2} + {\overset{\cdot}{u}}_{1}^{2}] {ds}_{1} + \frac{1}{2} \int_{0}^{l_{2}} ρ_{2} [{({\overset{\cdot}{v}}_{1} (l_{1}) + {\overset{\cdot}{u}}_{2} + \overset{\cdot}{w})}^{2} + {({\overset{\cdot}{u}}_{1} (l_{1}) - {\overset{\cdot}{v}}_{2})}^{2}] {ds}_{2} + \frac{1}{2} \int_{0}^{l_{3}} ρ_{3} {[{({\overset{\cdot}{v}}_{1} (l_{1}) + {\overset{\cdot}{u}}_{2} (l_{2}) - {\overset{\cdot}{v}}_{3} + \overset{\cdot}{w})}^{2} + ({\overset{\cdot}{u}}_{1} (l_{1}) - {\overset{\cdot}{v}}_{2} (l_{2}) - {\overset{\cdot}{u}}_{3})]}^{2} {ds}_{3} + \frac{1}{2} \int_{0}^{l_{4}} ρ_{4} [{({\overset{\cdot}{v}}_{1} (l_{1}) - {\overset{\cdot}{u}}_{4} + \overset{\cdot}{w})}^{2} + {({\overset{\cdot}{u}}_{1} (l_{1}) + {\overset{\cdot}{v}}_{4})}^{2}] {ds}_{4} \end{matrix}

(6)

where v_i is the transverse deflection, u_i is the axial deflection, and ρ_i is the linear mass density of the ith beam, and

\begin{matrix} U = \frac{1}{2} \int_{0}^{l_{1}} {EI}_{1} {(\frac{\partial^{2} v_{1}}{\partial s_{1}^{2}})}^{2} {ds}_{1} + \frac{1}{2} \int_{0}^{l_{2}} {EI}_{2} {(\frac{\partial^{2} v_{2}}{\partial s_{2}^{2}})}^{2} {ds}_{2} + \frac{1}{2} \int_{0}^{l_{3}} {EI}_{3} {(\frac{\partial^{2} v_{3}}{\partial s_{3}^{2}})}^{2} {ds}_{3} + \frac{1}{2} \int_{0}^{l_{4}} {EI}_{4} {(\frac{\partial^{2} v_{4}}{\partial s_{4}^{2}})}^{2} {ds}_{4} + \frac{1}{2} \int_{0}^{l_{1}} {EA}_{1} {(\frac{\partial u_{1}}{\partial s_{1}})}^{2} {ds}_{1} + \frac{1}{2} \int_{0}^{l_{2}} {EA}_{2} {(\frac{\partial u_{2}}{\partial s_{2}})}^{2} {ds}_{2} + \frac{1}{2} \int_{0}^{l_{3}} {EA}_{3} {(\frac{\partial u_{3}}{\partial s_{3}})}^{2} {ds}_{3} + \frac{1}{2} \int_{0}^{l_{4}} {EA}_{4} {(\frac{\partial u_{4}}{\partial s_{4}})}^{2} {ds}_{4} \end{matrix}

(7)

where E is the material's Youngs Modulus, I_i are the area moments of inertia, and A_i are the cross sectional areas of the beams.

It is assumed that the structure vibrates with finite amplitudes about its static equilibrium position. A small dimensionless parameter ε is introduced to help keep track of significant terms in the system response [7]. It is also assumed that the beams are inextensible and hence the inextensibility constraint can be written as [16]

{(1 + \frac{\partial u_{i}}{\partial s_{i}})}^{2} + {(\frac{\partial v_{i}}{\partial s_{i}})}^{2} = 1, i = 1, 2, 3

(8)

Assuming that the system response consists predominantly of the lowest two modes, the transverse deflection of the ith beam can be written as

v_{i} = ɛ (A_{1} φ_{i 1} + A_{2} φ_{i 2})

(9)

where A₁ and A₂ are the modal amplitudes that are assumed to be only functions of time, and ϕ_i₁ and ϕ_i₂ are the first and second mode shapes of the ith 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₁ and A₂ are nondimensional. Using the inextensibility constraint (Eq. ), the axial displacement can be obtained as [7]

u_{i} = ɛ^{2} (A_{1}^{2} η_{i 1} + A_{2}^{2} η_{i 2} + 2 A_{1} A_{2} η_{i 12}), i = 1, 2, 3, 4

(10)

where

η_{i 1} = - \frac{1}{2} \int_{0}^{s_{i}} {(\frac{d φ_{i 1}}{{ds}_{i}})}^{2} {ds}_{i}

(11a)

η_{i 2} = - \frac{1}{2} \int_{0}^{s_{i}} {(\frac{d φ_{i 2}}{{ds}_{i}})}^{2} {ds}_{i}

>

(11b)

η_{i 12} = - \frac{1}{2} \int_{0}^{s_{i}} (\frac{d φ_{i 1}}{{ds}_{i}}) (\frac{d φ_{i 2}}{{ds}_{i}}) {ds}_{i}

(11c)

The base excitation is also assumed to be of the form

\overset{··}{w} = ɛ^{2} F \cos (Ω t)

(12)

where Ω is the excitation frequency which is close to the second natural frequency of the structure. The difference between the excitation frequency and the second natural frequency is represented by the parameter σ₂ known as external mistuning

Ω = ω_{2} + ɛ σ_{2}

(13)

Similarly, another mistuning parameter, the internal mistuning (σ₁), is introduced to take into account the deviation of the lowest two natural frequencies from the perfect 1:2 ratio

ω_{2} = 2 ω_{1} + ɛ σ_{1}

(14)

To obtain the two-mode dynamic model of the system, and to further study the response using the method of averaging [12], the modal amplitudes are assumed to be of the form

A_{1} = p_{1} \cos (\frac{Ω}{2} t) + q_{1} \sin (\frac{Ω}{2} t)

(15a)

A_{2} = p_{2} \cos (Ω t) + q_{2} \sin (Ω t)

(15b)

where p_i and q_i are amplitude components which vary on a slow time scale τ = εt, as defined in Refs. [6,7].

The Lagrangian (L) of the structure can now be obtained by substituting Eqs. , , , and and their derivatives into the expressions of the kinetic and strain energies of the structure defined in Eqs. and . This Lagrangian is then averaged over the period of the oscillation

\frac{4 π}{Ω}

⁠. The amplitudes p_is and q_is 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 O(ε³) 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

〈 L 〉 = \int_{0}^{T} (T - U) dt

(16)

Substituting Eqs. and into the averaged Lagrangian (⁠

〈 L 〉

⁠), the following Euler–Lagrange equations can be obtained for the modal amplitudes p_is and q_is:

{p'}_{1} + ζ_{1} p_{1} + (\frac{σ_{1} + σ_{2}}{2}) q_{1} + 0.0169 ω_{1} (p_{2} q_{1} - p_{1} q_{2}) = 0

(17a)

{q'}_{1} + ζ_{1} q_{1} - (\frac{σ_{1} + σ_{2}}{2}) p_{1} + 0.0169 ω_{1} (p_{1} p_{2} + q_{1} q_{2}) = 0

(17b)

{p'}_{2} + ζ_{2} p_{2} + σ_{2} q_{2} - 2 \times 0.0068 ω_{2} p_{1} q_{1} = 0

(17c)

{q'}_{2} + ζ_{2} q_{2} - σ_{2} p_{2} + 0.0068 ω_{2} (p_{1}^{2} - q_{1}^{2}) = 9.5661 \times 10^{4} \frac{F}{ω_{2}}

(17d)

where 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 (ζ₁ and ζ₂). The modal damping coefficients (ζ₁ and ζ₂) were set equal to 0.1 arbitrarily as no experiments were conducted to estimate the real damping coefficients. Equations (17a)–(17d) 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 $a_{1} = \sqrt{p_{1}^{2} + q_{1}^{2}}$ ⁠) 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.

Fig. 11

View large Download slide

Illustration of the nonlinear response of structure 1

Fig. 12

View large Download slide

Amplitude of mode 1 excited through internal resonance for structures 1, 2, 3, and 4 when subjected to the same level of excitation. The solid line represents the stable branch whereas the dotted line represents the unstable branch.

Fig. 13

View large Download slide

Amplitude of mode 1 for all structures (structures 1, 2, 3, and 4) shown in Fig. 12 to clearly show a quantitative difference in response amplitudes for the four structures. Note that the structures 1 and 3 have nearly identical amplitudes. The solid line represents the stable branch whereas the dotted line represents the unstable branch.

Fig. 14

View large Download slide

Point of maximum deflection (marked with a black dot) for mode 1 for the four structures designed for 1:2 internal resonance (all dimensions are in meters)

Fig. 15

View large Download slide

Tip amplitudes for the four structures designed for 1:2 internal resonance

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.

References

1.

Nayfeh

,

A. H.

,

Younis

,

M. I.

, and

Abdel-Rahman

,

E. M.

,

2007

, “

Dynamic Pull-In Phenomenon in MEMS Resonators

,”

Nonlinear Dyn.

,

48

(

1–2

), pp.

153

–

163

.10.1007/s11071-006-9079-z

Google Scholar

Crossref

2.

Younis

,

M. I.

,

2011

,

MEMS Linear and Nonlinear Statics and Dynamics

(Microsystems), 1st ed.,

Springer

,

New York

.10.1007/978-1-4419-6020-7

3.

Rhoads

,

J. F.

,

Shaw

,

S.

, and

Turner

,

K.

,

2010

, “

Nonlinear Dynamics and Its Applications in Micro- and Nanoresonators

,”

J. Dyn. Syst., Meas., Control

,

132

(

3

), p.

034001

.10.1115/1.4001333

Google Scholar

Crossref

4.

Vyas

,

A.

,

Bajaj

,

A. K.

,

Raman

,

A.

, and

Peroulis

,

D.

,

2009

, “

A Microresonator Design Based on Nonlinear 1:2 Internal Resonance in Flexural Structural Modes

,”

J. Microelectromech.Sys.

,

18

(

3

), pp.

744

–

762

.10.1109/JMEMS.2009.2017081

Google Scholar

Crossref

5.

Vyas

,

A.

,

Peroulis

,

D.

, and

Bajaj

,

A. K.

,

2008

, “

Dynamics of a Nonlinear Microresonator Based on Resonantly Interacting Flexural-Torsional Modes

,”

Nonlinear Dyn.

,

54

(

1–2

), pp.

31

–

52

.10.1007/s11071-007-9326-y

Google Scholar

Crossref

6.

Bajaj

,

A. K.

,

Chang

,

S. I.

, and

Johnson

,

J. M.

,

1994

, “

Amplitude Modulated Dynamics of a Resonantly Excited Autoparametric Two Degree-of-Freedom System

,”

Nonlinear Dyn.

,

5

(

4

), pp.

433

–

457

.10.1007/BF00052453

Google Scholar

Crossref

7.

Balachandran

,

B.

, and

Nayfeh

,

A. H.

,

1990

, “

Nonlinear Motions of Beam-Mass Structure

,”

Nonlinear Dyn.

,

1

, pp.

39

–

61

.10.1007/BF01857584

Google Scholar

Crossref

8.

Wang

,

F.

, and

Bajaj

,

A.

,

2010

, “

Nonlinear Dynamics of a Three-Beam Structure With Attached Mass and Three-Mode Interactions

,”

Nonlinear Dyn.

,

62

(

1–2

), pp.

461

–

484

.10.1007/s11071-010-9734-2

Google Scholar

Crossref

9.

Saxena

,

A.

, and

Ananthasuresh

,

G. K.

,

2000

, “

On an Optimal Property of Compliant Topologies

,”

Struct. Multidiscip. Optim.

,

19

(

1

), pp.

36

–

49

.10.1007/s001580050084

Google Scholar

Crossref

10.

He

,

W.

,

Bindel

,

D.

, and

Govindjee

,

S.

,

2011

, “

Topology Optimization in Micromechanical Resonator Design

,”

Optim. Eng.

,

13

(2), pp.

271

–

292

.10.1007/s11081-011-9139-1

Google Scholar

Crossref

11.

Huang

,

X.

,

Zuo

,

Z. H.

, and

Xie

,

Y. M.

,

2010

, “

Evolutionary Topological Optimization of Vibrating Continuum Structures for Natural Frequencies

,”

Comput. Struct.

,

88

(

5–6

), pp.

357

–

364

.10.1016/j.compstruc.2009.11.011

Google Scholar

Crossref

12.

Nayfeh

,

A. H.

,

2000

,

Nonlinear Interactions: Analytical, Computational, and Experimental Methods

,

Wiley-Interscience

,

New York

.

13.

Rao

,

S. S.

,

2004

,

Mechanical Vibrations

, 4th ed.,

Prentice-Hall

,

Upper Saddle River, NJ

.

14.

The Ioffe Institute

,

2012

, “

Welcome to the Ioffe Institute

,” http://www.ioffe.ru/SVA/NSM/Semicond/Si/index.html

15.

Chong

,

E.

, and

Zak

,

S.

,

2008

,

An Introduction to Optimization

(Wiley Series in Discrete Mathematics and Optimization), 3rd ed.,

Wiley-Interscience

,

New York

.

16.

Nayfeh

,

A. H.

, and

Pai

,

P.

,

2000

,

Linear and Nonlinear Structural Mechanics

,

Wiley-Interscience

,

New York

.

Computational Synthesis for Nonlinear Dynamics Based Design of Planar Resonant Structures

Introduction

Linear Structure Synthesis

An Example.

Verification and Mode Shapes.

Results for Linear Synthesis.

Nonlinear Frequency Response Prediction

Discussion and Future Work

Acknowledgment

References

Get Email Alerts

Cited By

Journals

Conference Proceedings

eBooks

Resources

Opportunities

Computational Synthesis for Nonlinear Dynamics Based Design of Planar Resonant Structures

Introduction

Linear Structure Synthesis

An Example.

Verification and Mode Shapes.

Results for Linear Synthesis.

Nonlinear Frequency Response Prediction

Discussion and Future Work

Acknowledgment

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

Journals

Conference Proceedings

eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only