Dynamic analysis of a micro-resonator driven by electrostatic combs

Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Dynamic analysis of a micro-resonator driven by electrostatic combs M.T. Song a, D.Q. Cao a, W.D. Zhu a,b,⇑ a b

School of Astronautics, Harbin Institute of Technology, PO Box 137, Harbin 150001, China Department of Mechanical Engineering, University of Maryland Baltimore County, Baltimore, MD 21250, USA

a r t i c l e

i n f o

Article history: Received 9 September 2010 Received in revised form 4 December 2010 Accepted 5 December 2010 Available online 9 December 2010 Keywords: Micro-resonator Flexible beam Rigid body Global mode Orthogonality relation Repeated natural frequencies Gram–Schmidt orthogonalization method Nonlinear dynamic behavior

a b s t r a c t The dynamic response of a micro-resonator driven by electrostatic combs is investigated in this work. The micro-resonator is assumed to consist of eight flexible beams and three rigid bodies. The nonlinear partial differential equations that govern the motions of the flexible beams are obtained, as well as their boundary and matching conditions. The natural matching conditions for the flexible beams are the governing equations for the rigid bodies. The undamped natural frequencies and mode shapes of the linearized model of the microresonator are determined, and the orthogonality relation of the undamped global mode shapes is established. The modified Newton iterative method is used to simultaneously solve for the frequency equation and identify repeated natural frequencies that can occur in the micro-resonator and their multiplicities. The Gram–Schmidt orthogonalization method is extended to orthogonalize the mode shapes of the continuous system corresponding to the repeated natural frequencies. The undamped global mode shapes are used to spatially discretize the nonlinear partial differential equations of the micro-resonator. The simulation results show that the geometric nonlinearities of the flexible beams can have a significant effect on the dynamic response of the micro-resonator. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Comb drive actuators, which have been demonstrated to be more effective and have a higher quality factor than parallelplate-capacitor drives [1], have been widely used in micro-electro-mechanical systems [2–6]. Usually, the folded beams in the comb drive actuators are modeled as springs, whose stiffnesses can be calculated by empirical formulae [1,7–9]. However, some dynamic characteristics of micro-resonators driven by electrostatic combs cannot be captured by the spring-mass models, and there are several advantages to model the folded beams as a continuum. First, a continuous system model can more closely represent the physical system of a micro-resonator, and one can determine in a predictive manner the stiffness properties of a micro-resonator. Second, a micro-resonator has more than one resonant frequency in each direction, and a spring-mass model can only predict one resonant frequency in each direction. The existence of multiple natural frequencies can result in rich dynamic phenomena especially when one considers the nonlinear factors. For instance, there can be nonlinear interactions between different modes. Third, while a micro-resonator usually works around its first natural frequency, higher modes can contribute in the dynamic response of the micro-resonator, especially in control of the micro-resonator, which requires a sufficient number of higher modes to be included in a modal expansion. The higher modes can also be included in the modal expansion to examine the convergence of the solution. Fourth, one can include the thermal stress effect in the stress–strain relation in a continuous system model and take the temperature field into consideration. ⇑ Corresponding author at: Department of Mechanical Engineering, University of Maryland Baltimore County, Baltimore, MD 21250, USA. Tel.: +1 410 455 3394; fax: +1 410 455 1052. E-mail addresses: [email protected] (M.T. Song), [email protected] (D.Q. Cao), [email protected] (W.D. Zhu). 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.12.004

3426

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

Nomenclature A Af Ap b bf c E EI f Fe g h hs I l lf lff l r li ; li ‘li ; ‘ri mP mTL ; mTR f M M

cross-sectional area of flexible beams area of fingers in electrostatic combs for damping calculation area of the whole plate with movable fingers flexible beam width finger width damping coefficient Young’s modulus of flexible beams bending stiffness of flexible beams temporal separable solution electric field force gap between movable and stationary fingers thickness of the micro-resonator gap between the movable structure and the substrate identity matrix flexible beam length finger length finger overlap local coordinates of the left and right ends of the ith flexible beam non-dimensional local coordinates of the left and right ends of the ith flexible beam mass of the whole plate with movable fingers masses of trusses extended mass matrix generalized mass matrix used to define the orthogonality relation between coefficient vectors corresponding to two global mode shapes nf number of movable fingers N number of included modes p multiplicity of natural frequencies qj generalized coordinates t time uj linear combination of mode shapes from the current method to represent the jth mode shape from ABAQUS for the case of repeated natural frequencies vj jth mode shape from ABAQUS for the case of repeated natural frequencies Va ac driving voltage Vd dc bias voltage displacement of the ith flexible beam in the lateral direction wi ~i w non-dimensional displacement of the ith flexible beam wTL ; wP ; wTR displacements of the rigid bodies TL, P, TR in the lateral direction W electric field energy between movable and stationary fingers x, y global coordinates xi local coordinate for the ith flexible beam ~xi non-dimensional local coordinate for the ith flexible beam e0 permittivity of vacuum er relative permittivity of air kj integer in the modified Newton iterative formula to solve for the jth natural frequency l dynamic viscosity m kinematic viscosity q mass density of the flexible beams s non-dimensional time /i mode shape of the ith flexible beam in the micro-resonator jth mode shape of the ith flexible beam in the micro-resonator /ij u1, u2 any two differentiable functions with respect to ~xi U global mode shape vector of the micro-resonator Uj jth global mode shape vector of the micro-resonator x0 parameter for non-dimensionalizing the time x natural frequency of the micro-resonator xF Fourier frequency in spectrum plots xj jth natural frequency from the current method xABA jth natural frequency from ABAQUS j ~ x non-dimensional natural frequency of the micro-resonator

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

3427

~j x jth non-dimensional natural frequency of the resonator  ~ j  ; k ¼ 0; 1; . . . iterative sequence in the modified Newton iterative method to solve for the jth natural frequency x k

X e X ð_Þ ()0 ()(4) ()T kk

frequency of the ac voltage non-dimensional frequency of the ac voltage differentiation with respect to t or s differentiation with respect to xi or ~xi fourth derivative with respect to xi or ~xi Transpose of a matrix Euclidean norm of a vector

The finite element method has been commonly used to analyze the dynamic behaviors of micro-resonators [10–13]. However, the use of a finite element code to simulate complex multi-body systems, such as a micro-resonator, is prohibitively cumbersome, expensive, and time consuming. Hence it is crucial to develop an effective modeling technique that can accurately capture the dynamic characteristics of a micro-resonator. One of the difficulties in modeling a micro-resonator consisting of flexible beams and rigid bodies is to handle the coupling between the flexible beams and rigid bodies. If the local modes of the flexible beams are to be used in the spatial discretization, they need to satisfy all the boundary and matching conditions for the flexible beams, which can pose a considerable challenge. In this work, the undamped global mode shapes of the linearized model of a micro-resonator are derived for the first time and used in a Galerkin procedure to spatially discretize the nonlinear partial differential equations of the micro-resonator. Only the motion of the micro-resonator in the lateral direction is considered. It is shown that even a one-degree-of-freedom model of the micro-resonator can be sufficiently accurate, which demonstrates the effectiveness of the modeling technique presented here. Some interesting nonlinear behaviors, such as jump phenomena and superharmonic resonances, are observed in the micro-resonator. While repeated natural frequencies can commonly occur in two- or three-dimensional motions of symmetrical structures [14], it is interesting to find that they can also occur in the one-dimensional motion of the micro-resonator. A new methodology is developed to identify repeated natural frequencies and their multiplicities and to orthogonalize the mode shapes corresponding to the repeated natural frequencies. The results from the current method are in excellent agreement with those from the finite element method using the commercial software ABAQUS. 2. Equations of motion Fig. 1 shows a typical comb-driven micro-resonator. The movable structure suspended by folded beams at two anchors is actuated in the lateral direction by applying a dc bias voltage Vd and two ac driving voltages denoted by Va(t) and Va(t), which have the same amplitude and opposite directions. The electric field energy between the static and movable fingers changes as the micro-resonator vibrates and one can calculate the electric field force through this change. To analyze the dynamics of the movable plate (i.e., the center plate and the side plates in Fig. 1) in the actuated direction, one can model the micro-resonator in Fig. 1 as a system consisting of the flexible beams 1, 2, . . . , 8 and the rigid bodies TL, P, TR, as shown in Fig. 2. The beams join the rigid bodies at the junctions S1, S2, . . . , S12. There is an external excitation force Fe, which is the electric field force, in the y direction on the rigid body P in the xy plane. Hence one can only consider the motion of the micro-resonator in the y direction.

Fig. 1. Schematic of a typical comb-driven micro-resonator.

3428

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

Fig. 2. Schematic of the micro-resonator model.

Neglecting the rotary inertia and shear deformation effects and considering the effect of the axial force on the lateral motion yields the equation of motion of the ith flexible beam [15] ð4Þ

€i  EIwi þ qAw

" Z r 2 # 2 EA li @wi @ wi dxi ¼ 0; 2l lli @xi @x2i

i ¼ 1; 2; . . . ; 8

ð1Þ

3

1 where I ¼ 12 hb , and the nonlinear term results from the geometric nonlinearity of the flexible beam. Since the rigid body TL is subjected to the shear forces of the flexible beams 1, 2, 3, 4 at the junctions S1, S2, S3, S4, respectively, the dynamic equation of the rigid body TL is

000 000 000 € EIw000 1 ð0; tÞ  EIw2 ðl; tÞ  EIw3 ðl; tÞ  EIw4 ð0; tÞ ¼ mTL wTL

ð2Þ

Similarly, the dynamic equations of the rigid bodies P and TR are 000 000 000 _ € EIw000 1 ðl; tÞ þ EIw4 ðl; tÞ  EIw5 ðl; tÞ  EIw8 ðl; tÞ  c wP þ F e ¼ mP wP

ð3Þ

000 000 000 € EIw000 5 ð0; tÞ þ EIw6 ðl; tÞ þ EIw7 ðl; tÞ þ EIw8 ð0; tÞ ¼ mTR wTR

ð4Þ

and

respectively. The electric field energy between the movable and stationary fingers is

W¼

1 nf e0 er hðlff  wP Þ 1 nf e0 er hðlff þ wP Þ ðV d  V a sinðXtÞÞ2 þ ðV d þ V a sinðXtÞÞ2 2 g 2 g

ð5Þ

The electric field force in the y direction is

Fe ¼

@W 2nf he0 er ¼ V d V a sinðXtÞ @y g

ð6Þ

The Stokes’ damping model can be used to calculate the damping coefficient of the rigid body P [9,16]:

  sinhð2vhs Þ þ sinð2vhs Þ Af sinhð2vgÞ þ sinð2vgÞ c ¼ cp vhs 1 þ þ coshð2vhs Þ  cosð2vhs Þ Ap coshð2vgÞ  cosð2vgÞ

ð7Þ

 1=2 lA where cp ¼ hsp and v ¼ 2Xm . The boundary conditions of the flexible beams at A, B, C, D are

A: B:

w2 ð0; tÞ ¼ 0; w3 ð0; tÞ ¼ 0;

w02 ð0; tÞ ¼ 0 w03 ð0; tÞ ¼ 0

C:

w6 ð0; tÞ ¼ 0;

w06 ð0; tÞ ¼ 0

ð10Þ

D:

w7 ð0; tÞ ¼ 0;

w07 ð0; tÞ ¼ 0

ð11Þ

ð8Þ ð9Þ

The geometric matching conditions, which relate the displacements and slopes at the connecting points of the flexible beams, are S1, S2, S3, S4:

w1 ð0; tÞ ¼ w2 ðl; tÞ ¼ w3 ðl; tÞ ¼ w4 ð0; tÞ

ð12Þ

w01 ð0; tÞ ¼ w02 ðl; tÞ ¼ w03 ðl; tÞ ¼ w04 ð0; tÞ ¼ 0

ð13Þ

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

3429

S5, S6, S7, S8:

w1 ðl; tÞ ¼ w4 ðl; tÞ ¼ w5 ðl; tÞ ¼ w8 ðl; tÞ w01 ðl; tÞ ¼ w04 ðl; tÞ ¼ w05 ðl; tÞ ¼ w08 ðl; tÞ ¼ 0

ð14Þ ð15Þ

S9, S10, S11, S12:

w5 ð0; tÞ ¼ w6 ðl; tÞ ¼ w7 ðl; tÞ ¼ w8 ð0; tÞ

ð16Þ

w05 ð0; tÞ ¼ w06 ðl; tÞ ¼ w07 ðl; tÞ ¼ w08 ð0; tÞ ¼ 0

ð17Þ

The natural matching conditions, which relate the bending moments and shear forces at the connecting points of the flexible beams, are obtained from Eqs. (2)–(4), noting that w1 ð0; tÞ ¼ wTL ; w1 ðl; tÞ ¼ wP , and w5 ð0; tÞ ¼ wTR : 000 000 000 €  EIw000 1 ð0; tÞ  EIw2 ðl; tÞ  EIw3 ðl; tÞ  EIw4 ð0; tÞ ¼ mTL w1 ð0; tÞ 000 000 000 000 _ 1 ðl; tÞ þ F e ¼ mP w € 1 ðl; tÞ EIw1 ðl; tÞ þ EIw4 ðl; tÞ  EIw5 ðl; tÞ  EIw8 ðl; tÞ  cw

ð19Þ

000 000 000 € EIw000 5 ð0; tÞ þ EIw6 ðl; tÞ þ EIw7 ðl; tÞ þ EIw8 ð0; tÞ ¼ mTR w5 ð0; tÞ

ð20Þ

ð18Þ

With the introduction of the non-dimensional variables:

xi ; l

~xi ¼

where x0 ¼

~i ¼ w

qffiffiffiffiffiffiffiffiffi EI

qhbl4

¼

wi ; l

e ¼ s ¼ x0 t; X

ffi qffiffiffiffiffiffiffiffi 2 Eb 12ql4

X

x0

is a measure of the fundamental natural frequency of a flexible beam of the micro-resonator, Eqs.

(1) and (18)–(20) become

EI

x

2 3 0l

EA €~ ~ i; w ~ i Þw ~ 00i ¼ 0; ~ ð4Þ Ci ðw w i þ qAlwi  2x20 l

i ¼ 1; 2;    ; 8

~ 000 ~ 000 ~ 000 ~ 000 w 1 ð0; sÞ  w2 ð1; sÞ  w3 ð1; sÞ  w4 ð0; sÞ ¼

mTL € ~ 1 ð0; sÞ w

c

c ~ 000 c ~ 000 c ~ 000 c ~ 000 Fe c _ €~ 1 ð1; sÞ ~ 1 ð1; sÞ ¼ w w w1 ð1; sÞ þ w ð1; sÞ  w ð1; sÞ  w ð1; sÞ þ  mP 4 mP 5 mP 8 mP x20 l x0 mP

mP

~ 000 ~ 000 ~ 000 ~ 000 w 5 ð0; sÞ þ w6 ð1; sÞ þ w7 ð1; sÞ þ w8 ð0; sÞ ¼

mTR € ~ 5 ð0; sÞ w

c

ð21Þ

ð22Þ

ð23Þ

ð24Þ

where

c¼

EI

x20 l

; 3

Ci ðu1 ; u2 Þ ¼

Z ‘li

‘ri

@ u1 @ u2 d~xi ; @ ~xi @ ~xi

Fe ¼

2nf he0 er e sÞ V d V a sinð X g

The non-dimensional boundary and matching conditions can be similarly obtained from Eqs. (8)–(17). 3. Undamped natural frequencies and mode shapes of the linearized model The undamped, linearized equations of Eq. (21) are 4€ ~ ~ ð4Þ w i þ a wi ¼ 0;

where a ¼



qAl4 x20 EI

14

i ¼ 1; 2; . . . ; 8

ð25Þ

~ i ð~ . Substituting the separable solution w xi ; sÞ ¼ /i ð~ xi Þf ðsÞ into Eq. (25) yields

€f ðsÞ þ x ~ 2 f ð sÞ ¼ 0

ð26Þ

and ð4Þ

~ 2 /i ¼ 0 /i  a4 x

ð27Þ

Similarly, the non-dimensional boundary and matching conditions reduce to

/2 ð0Þ ¼ 0;

/02 ð0Þ ¼ 0

ð28Þ

B:

/3 ð0Þ ¼ 0;

ð29Þ

C:

/6 ð0Þ ¼ 0;

D:

/7 ð0Þ ¼ 0;

/03 ð0Þ ¼ 0 /06 ð0Þ ¼ 0 /07 ð0Þ ¼ 0

A:

ð30Þ ð31Þ

3430

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

S1, S2, S3, S4:

/1 ð0Þ ¼ /2 ð1Þ ¼ /3 ð1Þ ¼ /4 ð0Þ

ð32Þ

/01 ð0Þ ¼ /02 ð1Þ ¼ /03 ð1Þ ¼ /04 ð0Þ ¼ 0

ð33Þ

S5, S6, S7, S8:

/1 ð1Þ ¼ /4 ð1Þ ¼ /5 ð1Þ ¼ /8 ð1Þ

ð34Þ

/01 ð1Þ ¼ /04 ð1Þ ¼ /05 ð1Þ ¼ /08 ð1Þ ¼ 0

ð35Þ

S9, S10, S11, S12:

/5 ð0Þ ¼ /6 ð1Þ ¼ /7 ð1Þ ¼ /8 ð0Þ /05 ð0Þ ¼ /06 ð1Þ ¼ /07 ð1Þ ¼ /08 ð0Þ ¼ 0

ð36Þ ð37Þ

000 000 000  /000 1 ð0Þ  /2 ð1Þ  /3 ð1Þ  /4 ð0Þ ¼  000 000 000 /000 1 ð1Þ þ /4 ð1Þ  /5 ð1Þ  /8 ð1Þ ¼  000 000 000 /000 5 ð0Þ þ /6 ð1Þ þ /7 ð1Þ þ /8 ð0Þ ¼ 

mTR

c

mTL

mP

c

c

~ 2 /1 ð0Þ x

~ 2 /1 ð1Þ x

ð38Þ ð39Þ

~ 2 /5 ð0Þ x

ð40Þ

The solution of Eq. (27) is

/i ð~xÞ ¼ Ai cosðb~xi Þ þ Bi sinðb~xi Þ þ C i coshðb~xi Þ þ Di sinhðb~xi Þ

ð41Þ

1 2

~ , and Ai, Bi, Ci, and Di are arbitrary constants. Let where b ¼ ax

wi ¼ ½Ai

Bi

W ¼ ½wT1

wT2

Di T

Ci

ð42Þ

and

wT3

wT4

wT5

wT6

wT7

wT8 T

ð43Þ

The undamped global mode shapes of the micro-resonator are

U ¼ ½/1

/2

/8  ¼ wT B



ð44Þ

where

2 6 6 B¼6 6 4

3

B1 B2 0

7 7 7 7 5

0 .. . B8

in which

Bi ¼ ½cosðb~xi Þ

sinðb~xi Þ

coshðb~xi Þ

sinhðb~xi ÞT ;

i ¼ 1; 2; . . . ; 8

Substituting the boundary conditions in Eqs. (28)–(31) and the matching conditions in Eqs. (32)–(40) into Eq. (41), and using Eqs. (42) and (43), yields

~ ÞW ¼ 0 Gðx

ð45Þ

~ Þ 2 R3232 are given in Appendix A; all of the undefined entries in Gðx ~ Þ are zero. The frewhere entries of the matrix Gðx ~ ÞÞ ¼ 0; it has an infinite number of roots x ~ 1; x ~ 2 ; . . ., which are the nonquency equation can be determined from detðGðx dimensional natural frequencies of the micro-resonator. The numerical solution of the frequency equation can be sought ~ ÞÞ versus x, and subsequently find the small sections in which the curve as follows: first, one can draw the curve of detðGðx  ~ ÞÞ ¼ 0. Next, one can select initial points x ~ j  ðj ¼ 1; 2; . . .Þ in the small sections and use the intersect with the line detðGðx 0 modified Newton iterative formula [17]





~ j ~ j   kj x ¼x kþ1 k

 ~ j  ÞÞ detðGðx k d  ÞÞÞ ~ ðdetðGð x j ~ dx

ð46Þ

k

~ j ðj ¼ 1; 2; . . .Þ. If a natural frequency has a multiplicity p(p P 1), the iteration in where kj are positive integers, to solve for x Eq. (46) yields a quadratic convergence when kj = p and a linear convergence when kj – p [17]. When xj and its multiplicity p are determined, one can select 32  p independent equations from Eq. (45) to obtain p basic solutions that can span the linear solution space. Then one can select any p linearly independent solutions in the linear solution space to obtain p orthogonal solutions by letting them satisfy the orthogonality relation for the mode shapes when p > 1 as discussed in Section 4.2.

3431

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

One can then use the p orthogonal solutions W in Eq. (44) to obtain p orthogonal mode shapes. The undamped global mode shapes corresponding to distinct and repeated natural frequencies are orthogonal with respect to the extended mass matrix, as shown in Section 4. 4. Orthogonality of undamped global mode shapes Eq. (27) is written as

~ 2j M Uj ¼ L Uj x

ð47Þ

where the matrix M and the operator L are defined by

M ¼ qAlI;

"

@4 @4 @4 L¼ diag ; ; . . . ; @ ~x41 @ ~x42 @ ~x48 x20 l3 EI

# ð48Þ

4.1. The case of distinct natural frequencies (p = 1) Denote the two undamped global mode shapes corresponding to two distinct natural frequencies by Ur ð~ xÞ and Us ð~ xÞ, respectively. By Eq. (47)

~ 2r M Ur ¼ L Ur x

ð49Þ

Pre-multiplying Eq. (49) by UTs ð~ xÞ yields

~ 2r UTs M Ur ¼ UTs L Ur x

ð50Þ

Integrating Eq. (50) with respect to ~ xi from ‘li to ‘ri for the ith flexible beam (i = 1, 2, . . . , 8) yields 8 Z X i¼1

‘ri

‘li

~ 2r /is /ir d~xi ¼ qAlx

8 Z X i¼1

‘ri

EI

x

‘li

2 3 0l

/is

@ 4 /ir d~xi @ ~x4i

ð51Þ

Integration by parts and using the boundary conditions in Eqs. (28)–(31) and the matching conditions in Eqs. (32)–(40) yields 8 Z X i¼1

‘ri

‘li

~ 2r /is /ir d~xi ¼ qAlx

8 Z X i¼1

‘li

‘ri

EI

x

2 3 0l

/is

@ 4 /ir d~xi @ ~x4i

~ 2r mP /1s ð1Þ/1r ð1Þ  x ~ 2r mTL /1s ð0Þ/1r ð0Þ  x ~ 2r mTR /5s ð0Þ/5r ð0Þ þ ¼ x

Z 8 X EI i¼1

x

2 3 0l

‘ri

‘li

/00is /00ir d~xi

ð52Þ

/00is /00ir d~xi

ð53Þ

Exchanging the subscripts s and r in Eq. (52) yields 8 Z X i¼1

‘ri

‘li

~ 2s /is /ir d~xi ¼ x ~ 2s mTL /1s ð0Þ/1r ð0Þ  x ~ 2s mP /1s ð1Þ/1r ð1Þ  x ~ 2s mTR /5s ð0Þ/5r ð0Þ þ qAlx

Z 8 X EI i¼1

x20 l3

‘li

‘ri

Subtracting Eq. (53) from Eq. (52) yields



~ 2r

~ 2s

x x

" 8  X

qAl

i¼1

Z ‘li

‘ri

# /is /ir d~xi þ mTL /1s ð0Þ/1r ð0Þ þ mP /1s ð1Þ/1r ð1Þ þ mTR /5s ð0Þ/5r ð0Þ ¼ 0

ð54Þ

Let

e rs ¼ M

8 X i¼1

qAl

Z ‘li

‘ri

/is /ir d~xi þ mTL /1s ð0Þ/1r ð0Þ þ mP /1s ð1Þ/1r ð1Þ þ mTR /5s ð0Þ/5r ð0Þ

ð55Þ

~r – x ~ s and by Eqs. (54) and (55) be entries of the extended mass matrix f M. When r – s; x

e rs ¼ 0 M

ð56Þ

Hence Ur and Us are orthogonal with respect to M. 4.2. The case of repeated natural frequencies (p > 1) ~k ¼ x ~ kþ1 ¼    ¼ x ~ kþp1 , one can extend the If there exist repeated natural frequencies with a multiplicity p, i.e., x Gram–Schmidt orthogonalization method for a discrete system [18] to orthogonalize the corresponding mode shapes for the continuous system here. To this end, Eq. (56) is written as

3432

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

 Ws ¼ 0 WTr M

ð57Þ

where entries of the generalized mass matrix M are given in Appendix B. Note that Eq. (57) for the continuous system here has a form similar to that for a discrete system. Let Wk, Wk+1, . . . , Wk+p1 be a set of linearly independent solutions of Eq. (45) ~k ¼ x ~ kþ1 ¼    ¼ x ~ kþp1 . Define a set of new vectors Wk ; Wkþ1 ; . . . ; Wkþp1 : corresponding to the repeated natural frequencies x

8 W k ¼ Wk > > > > > W ¼ Wkþ1 þ d1;1 Wk > > < kþ1 Wkþ2 ¼ Wkþ2 þ d2;1 Wk þ d2;2 Wkþ1 > > .. > > > . > > : Wkþp1 ¼ Wkþp1 þ dp1;1 Wk þ dp1;2 Wkþ1 þ    þ dp1;p1 Wkþp2

ð58Þ

where the coefficients d1,1, . . . , dp1,p1 are determined in order to satisfy the orthogonality relation

WTr MWs ¼ 0 ðr – s; r; s ¼ k; k þ 1;    ; k þ p  1Þ

ð59Þ

WTk MWkþ1 ¼ WTk M Wkþ1 þ d1;1 WTk MWk ¼ 0

ð60Þ

From

one has

d1;1 ¼ 

WTk M Wkþ1

ð61Þ

WTk MWk

From

(

WTk MWkþ2 ¼ WTk M Wkþ2 þ d2;1 WTk MWk ¼ 0

ð62Þ

WTkþ1 MWkþ2 ¼ WTkþ1 M Wkþ2 þ d2;2 WTkþ1 MWkþ1 ¼ 0 one has

d2;1 ¼ 

WTk M Wkþ2 WTk MWk

;

d2;2 ¼ 

WTkþ1 M Wkþ2

ð63Þ

WTkþ1 MWkþ1

From

8 T Wk MWkþp1 ¼ WTk M Wkþp1 þ dp1;1 WTk MWk ¼ 0 > > > > > < WTkþ1 MWkþp1 ¼ WTkþ1 M Wkþp1 þ dp1;2 WTkþ1 MWkþ1 ¼ 0 .. > > . > > > : T Wkþp2 MWkþp1 ¼ WTkþp2 M Wkþp1 þ dp1;p1 WTkþp2 MWkþp2 ¼ 0

ð64Þ

one has

dp1;1 ¼ 

WTk M Wkþp1 WTk MWk

;

dp1;2 ¼ 

WTkþ1 M Wkþp1 WTkþ1 MWkþ1

;

...

;

dp1;p1 ¼ 

WTkþp2 M Wkþp1 WTkþp2 MWkþp2

ð65Þ

Note that the orthogonal vectors Wk ; Wkþ1 ; . . . ; Wkþp1 cannot be uniquely determined. The forms of Wk ; Wkþ1 ; . . . ; Wkþp1 depend on the linearly independent solutions Wk, Wk+1, . . . , Wk+p1 selected from the p-dimensional linear solution space of Eq. (45), which are not unique, and there are an infinite number of ways to subsequently find a set of orthogonal mode shapes. 5. Spatial discretization The displacements of the flexible beams are approximated by a finite series using the undamped global mode shapes determined in Sections 3 and 4.2:

~ i ð~xi ; sÞ ¼ w

N X

/ij ð~xi Þqj ðsÞ

ð66Þ

j¼1

Substituting Eq. (66) into Eq. (21) yields N EI X 2 3 0 l j¼1

x

ð4Þ /ij ð~xÞqj ðsÞ þ qAl

N X j¼1

€ j ð sÞ  /ij ð~xÞq

N X N X N EA X Ci ð/ij ; /ik Þ/00im qj qk qm ¼ 0 2 2x0 l m¼1 k¼1 j¼1

ð67Þ

3433

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

Using Eq. (27) in Eq. (67) yields

qAl

N X

~ 2j /ij ð~xÞqj ðsÞ þ qAl x

N X

j¼1

€j ðsÞ  /ij ð~xÞq

j¼1

N X N X N EA X Ci ð/ij ; /ik Þ/00im qj qk qm ¼ 0 2x20 l m¼1 k¼1 j¼1

ð68Þ

Substituting Eq. (66) into Eqs. (22)–(24) yields N mTL X

c

j¼1

N mP X

c

€j ðsÞ þ /1j ð0Þq

/000 1j ð0Þqj ðsÞ þ

j¼1

€j ðsÞ þ /1j ð1Þq

j¼1

N X

/000 2j ð1Þqj ðsÞ þ

j¼1

N c X

cx0

N X

/1j ð1Þq_ j ðsÞ 

j¼1

N X

/000 3j ð1Þqj ðsÞ þ

N X

j¼1

N X

/000 1j ð1Þqj ðsÞ 

j¼1

N X

/000 4j ð0Þqj ðsÞ ¼ 0

ð69Þ

j¼1

/000 4j ð1Þqj ðsÞ þ

j¼1

N X

/000 5j ð1Þqj ðsÞ þ

j¼1

N X

/000 8j ð1Þqj ðsÞ 

j¼1

Fe

cx20 l

¼0 ð70Þ

N mTR X

c

€j ðsÞ  /5j ð0Þq

j¼1

N X

/000 5j ð0Þqj ðsÞ 

j¼1

N X

/000 6j ð1Þqj ðsÞ 

j¼1

N X

/000 7j ð1Þqj ðsÞ 

j¼1

N X

/000 8j ð0Þqj ðsÞ ¼ 0

ð71Þ

j¼1

Multiplying Eq. (68) by /in ð~ xi Þ and integrating the equation with respect to ~ xi from ‘li to ‘ri , multiplying Eq. (69) by /1n(0), multiplying Eq. (70) by /1n(1), multiplying Eq. (71) by /5n (0), adding all of the resulting equations, and using the orthogonality relations in Eqs. (56) and (59) yields

Z ‘r N 8 X N X N X N X i c X EA e nn q €n ðsÞ þ M /1j ð1Þ/1n ð1Þq_ j ðsÞ þ K nn qn ðsÞ  C ð/ ; / Þq q q /00im ð~xÞ/in ð~xÞd~x i ij ik j k m 2 l x0 j¼1 x 2l ‘ 0 i¼1 m¼1 k¼1 j¼1 i ¼

2nf he0 er ~ sÞ/ ð1Þ; V d V a sinðX 1n glx20

n ¼ 1; 2; . . .

ð72Þ

e nn . Note that while the mass and stiffness matrices in Eq. (72) have been diagonalized due to orthogonality ~ 2n M where K nn ¼ x of the undamped global mode shapes, the damping and nonlinear terms remain coupled, as expected. Eq. (72) can be numerically integrated using the Runge–Kutta method. 6. Results and discussion The values of the parameters of the micro-resonator in Ref. [7], as shown in Table 1, are used. The first several undamped natural frequencies and mode shapes of the micro-resonator are calculated and shown in Table 2. It is interesting to find that repeated natural frequencies occur due to symmetry of the micro-resonator. The fourth through eighth natural frequencies in Table 2 are repeated natural frequencies with a multiplicity of five, and the mode shapes shown there are the orthogonalized ones. While higher natural frequencies are not shown here, repeated natural frequencies with the same multiplicity occur after every three distinct natural frequencies. The undamped natural frequencies and mode shapes are also calculated using ABAQUS for validation purposes, as shown in Table 3. Each flexible beam is divided into 151 elements. The maximum relative error between the natural frequencies

Table 1 Values of the parameters of the micro-resonator. Parameter

Value

Finger gap g Finger length lf Finger overlap lff Finger width bf Movable finger numbers nf Flexible beam length l Flexible beam width b Dimensions of the center plate Dimensions of the side plate 1 2 Dimensions of the side plate 2 2 Thickness h Substrate gap hs Truss length Truss width Young’s modulus E Mass density q

2.88 lm 40.05 lm 19.44 lm 3.27 lm 30 151 lm 1.1 lm 54.9  19.26 lm2 89.6  14.13 lm2 175.47  11.3 lm2 1.96 lm 2 lm 78 lm 13 lm 150 GPa 2300 kg/m3

3434

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442 Table 2 Undamped natural frequencies and mode shapes of the micro-resonator. Natural frequency (rad/s)

Mode shape

x1 = 63,551

x2 = 282,651

x3 = 291,646

x4 = 2,516,248

x5 = 2,516,248

x6 = 2,516,248

x7 = 2,516,248

x8 = 2,516,248

x9 = 2,522,177

x10 = 2,639,790

from the current and finite element methods, defined by

    xj xABA  j xj

, is 0.034%. Due to non-uniqueness of the mode shapes cor-

responding to the repeated natural frequencies, the fourth through eighth mode shapes in Table 2 look different from those in Table 3, but correspondence between these mode shapes can be established. One can first discretize the fourth through eighth mode shapes from the current method using the same mesh as the finite element method. One approach is then to represent the fourth through eighth mode shapes from the finite element method, denoted by v4, v5, . . . , v8, as a linear combination of the corresponding discretized mode shapes from the current method, with the coefficients in the linear combination determined using the least squares method. The resulting linear combination of the discretized mode shapes from the current method, denoted by u4, u5, . . . , u8, are in excellent agreement with the corresponding mode shapes from the finite element method (v4, v5, . . ., v8), as shown in Table 4. The errors defined by

kv j uj k kv j k

are 0.0401%, 0.0083%, 0.6126%, 0.5672%,

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

3435

Table 3 Undamped natural frequencies and mode shapes of the micro-resonator from ABAQUS. Natural frequency (rad/s)

Mode shape

xABA ¼ 63; 547 1

xABA ¼ 282; 629 2

xABA ¼ 291; 624 3

xABA ¼ 2; 515; 426 4

xABA ¼ 2; 515; 426 5

xABA ¼ 2; 515; 426 6

xABA ¼ 2; 515; 426 7

xABA ¼ 2; 515; 426 8

xABA ¼ 2; 521; 319 9

xABA 10 ¼ 2; 638; 897

0.4900% for j = 4, 5, . . . , 8, respectively. The modal assurance criterion numbers [19] between the mode shapes from the current and finite element methods, as shown in Tables 2 and 3 for the first, second, third, ninth, and tenth modes and in Table 4 for the fourth through eighth modes, are all essentially 100%. One can see from Table 2 that the second natural frequency is much larger than the first natural frequency. Hence one can retain only the first mode (i.e., N = 1 in Eq. (66)) in calculating the dynamic response of the micro-resonator driven by the electric field force when X varies near x1. In addition, the fourth natural frequency is much larger than the third natural frequency, and one can retain the first three modes in Eq. (66). The results from the one-degree-of-freedom model will be compared with those from the three-degree-of-freedom model. In the following calculations, Vd is assumed to be 40 V,

3436

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442 Table 4 The linear combination of the fourth through eighth mode shapes determined using the current method: u4, u5, . . . ,u8. Natural frequency

Mode shape

x4

x5

x6

x7

x8

and the transient responses of the rigid body P are calculated for two thousand periods of the electric field force to ensure that the steady-state responses are reached. When Va varies from 0.1 to 50 V and X varies from 1,000 to 150,000 rad/s, the amplitude of the steady-state response of the rigid body P is shown in Fig. 3; the initial conditions are assumed to be zero, and the results obtained correspond to the case with increasing Va and decreasing X. The frequency response curves for Va = 1, 10, 20 V are shown in Fig. 4. The cases with increasing and decreasing X are obtained as follows: the initial conditions for the first point (X = 1,000 rad/s for the increasing X case and X = 250,000 rad/s for the decreasing X case) are assumed to be zero, and some steady-state response data at a current X are used as the initial conditions for the next X. It is seen that due to the geometric nonlinearities, the peaks of these curves do not exactly occur at the first undamped natural frequency of the linearized model of the micro-resonator and the peak frequencies depend on Va. A larger Va results in a larger peak frequency and a larger response amplitude. When Va approaches zero, the peak frequency approaches the first undamped natural frequency of the linearized model of

Fig. 3. The amplitude of the steady-state response of the rigid body P as a function of the amplitude and frequency of the ac voltage; N = 1.

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

3437

Fig. 4. The frequency response curves of the rigid body P for different amplitudes of the ac voltage; N = 1. An expanded view of the superharmonic resonance corresponding to X  15 x1 when va = 20 V is also shown.

the micro-resonator. The jumps for the cases with increasing and decreasing X occur at different X. In addition, superharmonic resonance corresponding to X  13 x1 is observed when Va = 10, 20 V, and that corresponding to X  15 x1 can also be observed when Va = 20 V. The geometric nonlinearity effect decreases as Va decreases. The frequency response curve from the one-degree-of-freedom model is in excellent agreement with that from the three-degree-of-freedom model, as shown in Fig. 5. When Va varies from 0.1 to 50 V, the force-response curves for X = 50,000, 70,000, 90,000 rad/s are shown in Fig. 6. The cases with increasing and decreasing Va are obtained as follows: the initial conditions for the first point (Va = 0.1 V for the increasing Va case and Va = 50 V for the decreasing Va case) are assumed to be zero, and some steady-state response data at a current Va are used as the initial conditions for the next Va. The jump phenomenon is more apparent for a larger X, and the jumps for the cases with increasing and decreasing Va occur at different Va. The force-response curve from the one-degree-of-freedom model is in excellent agreement with that from the three-degree-of-freedom model, as shown in Fig. 7. When X = 90,000 rad/s and Va = 25 V, the steady-state response of the rigid body P, calculated from the one-degree-offreedom model is shown in Fig. 8; the initial conditions are assumed to be zero. From the Poincaré map in Fig. 8(c), which is constructed by drawing a dot at all integer multiples of the period of the ac voltage, one can conclude that the rigid body P has a periodic steady-state motion. The steady-state response of the rigid body P, the phase portrait, the Poincaré map, and

Fig. 5. The frequency response curves of the rigid body P from the one-and three-degree-of-freedom models; Va = 20 V. An expanded view of the superharmonic resonance corresponding to X  15 x1 is also shown.

3438

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

Fig. 6. The force-response curves of the rigid body P for various values of X; N = 1.

Fig. 7. The force-response curves of the rigid body P from the one- and three-degree-of-freedom models; X = 90,000 rad/s.

the spectrum plot of the steady-state response from the three-degree-of-freedom model with the same initial conditions as those for the one-degree-of-freedom model, as shown in Fig. 9, are in excellent agreement with those from the one-degreeof-freedom model in Fig. 8. The small peaks that occur at 270,000 rad/s in Figs. 8(d) and 9(d) are associated with the thirdorder harmonics due to the cubic nonlinearities. 7. Conclusions A micro-resonator consisting of eight flexible beams and three rigid bodies is governed by eight nonlinear partial differential equations for the flexible beams with geometric and natural matching conditions, with the natural matching conditions being the governing equations for the rigid bodies. The undamped natural frequencies and mode shapes of the micro-resonator are exactly derived from the linearized model. An extended mass matrix is defined to account for the masses of the rigid bodies in the orthogonality relation for the mode shapes. Repeated natural frequencies with a multiplicity of five occur after every three distinct natural frequencies. The modified Newton iterative method can be effectively used to simultaneously solve for the frequency equation and identify the repeated natural frequencies and their multiplicities by examining its convergence rates. The Gram–Schmidt orthogonalization method for a discrete system is extended to the continuous system here to orthogonalize the mode shapes corresponding to the repeated natural frequencies. The calculated natural frequencies and mode shapes using the current method are in excellent agreement with those from the finite

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

3439

Fig. 8. (a) The steady-state response of the rigid body P, (b) the phase portrait, (c) the Poincaré map, and (d) the spectrum plot of the steady-state response; Va = 25 V, Vd = 40 V, X = 90,000 rad/s, and N = 1. An expanded view of the small peak in the spectrum in (d) is also shown.

Fig. 9. (a) The steady-state response of the rigid body P, (b) the phase portrait, (c) the Poincaré map, and (d) the spectrum plot of the steady-state response; Va = 25 V, Vd = 40 V, X = 90,000 rad/s, and N = 3. An expanded view of the small peak in the spectrum in (d) is also shown.

element method using ABAQUS. Correspondence between the mode shapes corresponding to the repeated natural frequencies from the current and finite element methods is established using the least squares method. The undamped, global mode shapes of the micro-resonator are used to spatially discretize the nonlinear partial differential equations. Use of the one- and three-degree-of-freedom models for the micro-resonator results in accurate response solutions. The geometric nonlinearities of the flexible beams can have a significant effect on the dynamic response of the micro-resonator. The methodology can be used to analyze similar systems in other applications.

Acknowledgements This work is supported by the National Natural Science Foundation of China under Grant No. 10772056, the ChangJiang Scholar Program of the Ministry of Education of China, the Harbin Institute of Technology, and the National Science

3440

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

Foundation through Award CMMI-0348605. The authors would like to thank Prof. Baisheng Wu for suggesting the use of the modified Newton iterative method to identify repeated natural frequencies. Appendix A. Entries of the matrix G in Eq. (45) Entries of the matrix G in Eq. (45) are:

g 1;5 ¼ 1;

g 1;7 ¼ 1

g 2;6 ¼ 1;

g 2;8 ¼ 1

g 3;9 ¼ 1;

g 3;11 ¼ 1

g 4;10 ¼ 1;

g 4;12 ¼ 1

g 5;21 ¼ 1;

g 5;23 ¼ 1

g 6;22 ¼ 1;

g 6;24 ¼ 1

g 7;25 ¼ 1;

g 7;27 ¼ 1

g 8;26 ¼ 1;

g 8;28 ¼ 1

g 9;1 ¼ x2 ;

g 9;2 ¼

g 9;8 ¼  g 9;14 ¼

cb3 cosh b

cb3 mTL

mTL

cb3 ; mTL

g 9;9 ¼

;

g 9;16 ¼ 

;

g 10;2 ¼ 1;

g 10;4 ¼ 1

g 11;1 ¼ 1;

g 11;3 ¼ 1;

g 12;5 ¼ sin b; g 13;1 ¼ 1;

g 13;3 ¼ 1;

g 16;14 ¼ 1;

mTL

cb3 ; mTL

g 9;10 ¼

;

g 9;5 ¼

cb3 cos b mTL

;

cb3 sin b mTL

g 9;11 ¼

g 9;6 ¼

;

cb3 sinh b mTL

;

cb3 cos b mTL

g 9;7 ¼

;

g 9;12 ¼ 

cb3 sinh b

cb3 cosh b mTL

mTL ;

cb3 mTL

g 11;5 ¼  cos b;

g 11;6 ¼ sin b;

g 12;7 ¼  sinh b;

g 13;9 ¼  cos b;

g 14;10 ¼ cos b;

g 15;3 ¼ 1;

g 9;4 ¼ 

cb3 sin b

g 12;6 ¼ cos b;

g 14;9 ¼ sin b; g 15;1 ¼ 1;

g 9;3 ¼ x2 ;

g 13;10 ¼ sin b;

g 14;11 ¼  sinh b;

g 15;13 ¼ 1;

g 11;7 ¼  cosh b;

g 11;8 ¼ sinh b

g 12;8 ¼ cosh b g 13;11 ¼  cosh b;

g 13;12 ¼ sinh b

g 14;12 ¼ cosh b

g 15;15 ¼ 1

g 16;16 ¼ 1

cb3 sin b

þ x2 cos b;

cb3 cos b

g 18;1 ¼ b sin b; g 19;1 ¼ cos b;

g 20;13 ¼ b sin b; g 21;1 ¼ cos b;

g 18;2 ¼ b cos b;

g 19;2 ¼ sin b;

g 19;13 ¼  cos b;

g 17;2 ¼ 

g 19;3 ¼ cosh b;

g 20;14 ¼ b cos b; g 21;18 ¼ sin b;

g 22;17 ¼ b sin b;

g 22;18 ¼ b cos b; g 23;2 ¼ sin b;

g 23;30 ¼ sin b;

g 24;29 ¼ b sin b;

g 24;30 ¼ b cos b;

g 18;4 ¼ b cosh b

g 19;4 ¼ sinh b g 19;16 ¼  sinh b

g 20;15 ¼ b sinh b;

g 20;16 ¼ b cosh b

g 21;4 ¼ sinh b

g 21;19 ¼  cosh b; g 22;19 ¼ b sinh b;

g 23;3 ¼ cosh b;

g 23;29 ¼  cos b;

g 17;3 ¼

g 19;15 ¼  cosh b;

g 21;3 ¼ cosh b;

g 21;17 ¼  cos b; g 23;1 ¼ cos b;

g 18;3 ¼ b sinh b;

g 19;14 ¼  sin b;

g 21;2 ¼ sin b;

þ x2 sin b;

cb3 sinh b

þ x2 cosh b; mP mP mP cb3 cosh b cb3 sin b cb3 cos b cb3 sinh b þ x2 sinh b; g 17;13 ¼ ; g 17;14 ¼  ; g 17;15 ¼ ; g 17;4 ¼ mP mP mP mP 3 3 3 3 cb cosh b cb sin b cb cos b cb sinh b cb3 cosh b ; g 17;17 ¼ ; g 17;18 ¼ ; g 17;19 ¼ ; g 17;20 ¼  ; g 17;16 ¼ mP mP mP mP mP cb3 sin b cb3 cos b cb3 sinh b cb3 cosh b ; g 17;30 ¼ ; g 17;31 ¼ ; g 17;32 ¼  g 17;29 ¼ mP mP mP mP

g 17;1 ¼

g 21;20 ¼ sinh b g 22;20 ¼ b cosh b

g 23;4 ¼ sinh b

g 23;31 ¼  cosh b; g 24;31 ¼ b sinh b;

g 23;32 ¼ sinh b g 24;32 ¼ b cosh b

;

3441

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

g 25;17 ¼ x2 ; g 25;22 ¼ 

g 25;18 ¼ 

cb3 cos b mTR

cb3 mTR

g 25;19 ¼ x2 ;

;

g 25;23 ¼

;

3

g 25;26 ¼ 

cb cos b mTR

g 25;27 ¼

;

g 26;20 ¼ 1

g 27;17 ¼ 1;

g 27;19 ¼ 1;

g 29;17 ¼ 1;

mTR

g 25;24 ¼

cb sinh b

g 28;22 ¼ cos b;

g 31;17 ¼ 1;

g 31;19 ¼ 1; g 32;32 ¼ 1

;

g 25;28 ¼

g 25;21 ¼

cb3 cosh b mTR

cb cosh b mTR

;

cb3 sin b mTR

g 25;25 ¼

;

cb3 sin b mTR

;

3

;

g 25;30 ¼ 

cb

mTR

g 27;23 ¼  cosh b;

;

g 25;32 ¼

cb3 mTR

g 27;24 ¼  sinh b

g 28;24 ¼ cosh b

g 29;26 ¼  sin b;

g 30;27 ¼ sinh b;

g 31;29 ¼ 1;

mTR

;

g 27;22 ¼  sin b;

g 28;23 ¼ sinh b;

g 29;25 ¼  cos b;

g 30;26 ¼ cos b;

g 32;30 ¼ 1;

mTR

cb3

3

g 27;21 ¼  cos b;

g 29;19 ¼ 1;

g 30;25 ¼  sin b;

;

3

g 26;18 ¼ 1;

g 28;21 ¼  sin b;

cb3 sinh b

g 25;20 ¼

g 29;27 ¼  cosh b;

g 29;28 ¼  sinh b

g 30;28 ¼ cosh b

g 31;31 ¼ 1

Appendix B. Entries of the matrix M in Eq. (57) The generalized mass matrix M in Eq. (57) is expressed as

2 6 6 6 6 6 6 6 M¼6 6 6 6 6 6 6 4

3

A1 þ B þ C

7 7 7 7 7 7 7 7 7 7 7 7 7 7 5

A2 A3

0 A4 A5 þ D

0

A6 A7 A8

where

3 cosðbj ~xi Þ 6 sinðb x~i Þ 7 j 7 6 7½cosðbj ~xi Þ 6 4 coshðbj ~xi Þ 5 sinhðbj ~xi Þ 2

Ai ¼ qAl

2

Z

‘ri ‘li

mTL

6 0 6 B¼6 4 mTL 0

0 mTL

0

07 7 7 05

0

0

0 0

coshðbj ~xi Þ

sinhðbj ~xi Þd~xi

3

0 mTL

0

sinðbj ~xi Þ

3 cos bj 6 sin b 7 j 7 6 C ¼ mP 6 7 cos bj 4 cosh bj 5 2

sin bj

cosh bj

sinh bj



sinh bj 2

mTR

6 0 6 D¼6 4 mTR 0

0 mTR

0

3

0 mTR

07 7 7 05

0

0

0

0 0

References [1] Tang WC, Nguyen CT-C, Howe RT. Laterally driven polysilicon resonant microstructures. Sensors Actuators A 1989;20:25–32. [2] Tang WC. Electrostatic comb drive for resonant sensor and actuator application PhD dissertation Department of Electrical Engineering and Computer Sciences University of California, Berkeley, CA, USA; 1990.

3442

M.T. Song et al. / Commun Nonlinear Sci Numer Simulat 16 (2011) 3425–3442

[3] Boser BE, Howe RT. Surface micromachined accelerometers. IEEE J Solid-State Circuit 1996;31:366–75. [4] Clark WA, Howe RT, Horowitz R. Surface micromachined Z-axis vibratory rate gyroscope, Tech. Dig. IEEE Solid-State Sensor Actuator Workshop. Hilton Head Island, SC; 1996. pp 283–7. [5] Lee AP. Impact actuation of micromechanical structures PhD dissertation University of California Berkeley, CA, USA; 1992. [6] Lin LW, Howe RT, Pisano AP. Microelectromechanical filters for signal processing. J Microelectromech Syst 1998;7:286–94. [7] Ye WJ, Wang X, Hemmert W, Freeman D, White J. Air damping in laterally oscillating microresonators: a numerical and experimental study. J Microelectromech Syst 2003;12:557–66. [8] Huang W, Lu GY. Analysis of lateral instability of in-plane comb drive MEMS actuators based on a two-dimensional model. Sensors Actuators A 2004;113:78–85. [9] Lee KB, Pisano AP, Lin LW. Nonlinear behaviors of a comb drive actuator under electrically induced tensile and compressive stresses. J Micromech Microeng 2007;17:557–66. [10] Shi F, Ramesh P, Mukherjee S. Simulation methods for micro-electro- mechanical structures (MEMS) with application to a microtweezer. Comput Struct 1995;56:769–83. [11] Ostergaard D, Gyimesi M. Finite element based reduced order modeling of micro electro mechanical systems (MEMS) MSM; 2000. pp. 684-7. [12] Ostergaard D. Tooling up for microelectromechanical systems (MEMS). ANSYS Solutions 2000;2:22–6. [13] Hung ES, Yang Y-J, Senturia SD. Low-order models for fast dynamical simulation of MEMS microstructures. In: Int Conf on Solid-State Sensors and Actuators, Transducer’97; 1997. pp. 1101–4 [14] Mohan SJ, Pratap R. A natural classification of vibration modes of polygonal ducts based on group theoretic analysis. J Sound Vibration 2004;269:745–64. [15] Younis MI, Abdel-Rahman EM, Nayfeh A. A reduced-order model for electrically actuated microbeam-based MEMS. J Microelectromech Syst 2003;12:672–80. [16] Li G, Hughes H. Review of viscous damping in micro-machined structures. Micromachined Devices Components VI, Proc SPIE 2000;4176:30–46. [17] Gerald CF, Wheatley PO. Applied numerical analysis. 7th ed. Pearson/Addison-Wesley; 2004. [18] Ginsberg JH. Mechanical and structural vibrations: theory and applications. New York: Wiley; 2001. [19] Allemang RJ, Brown DL. A correlation coefficient for modal vector analysis. In: Proc of the 1st Int. Modal Analysis Conf, Orlando, Florida; 1982.