Modelling the electrostatic actuation of MEMS: state of the art 2005.

Modelling the electrostatic actuation of MEMS: state of the art 2005. A. Fargas Marquès, R. Costa Castelló and A.M. Shkel IOC-DT-P-2005-18 Setembre 2005 MODELING THE ELECTROSTATIC ACTUATION OF MEMS. STATE OF THE ART 2005 A. Fargas Marqu`es∗, R. Costa Castell´o∗, and A. M. Shkel+ ∗Institut d’Organitzaci´o i Control de Sistemes Industrials, IOC-UPC and +University of California at Irvine July 2005 Abstract Most of MEMS devices are actuated using electrostatic forces. Parallel or lateral plate actuators are the types commonly used. Nevertheless, electrostatic actuation has some limitations due to its non-linear nature. This work presents a methodic overview of the existing techniques applied to the Micro-Electro-Mechanical Systems (MEMS) electrostatic actuation modeling and their implications to the dynamic behavior of the electromechanical system. 1 Introduction The field of Micro-Electro-Mechanical Systems (MEMS) has undergone a startling revolution in recent years. It is now possible to produce accelerometers less than one millimeter on a side, functioning motors that can only be seen with the aid of a microscope, gears smaller than a human hair, and needles so tiny they can deliver an injection without stimulating nerve cells. The use of existing integrated circuit technology in the design and production of MEMS devices allows these devices to be batch-manufactured, what in turn converts them due to their quantity in almost inexpensive. The first sector to benefit from this revolution has been the automotive industry, where devices and applications that once could only be dreamed about have suddenly been made possible and are used everywhere. The ability to manufacture mechanical parts such as resonators, sensors, gears and levers on a micron length scale is not however the end of the story. The challenge is also to understand and control the physical systems behavior on these scales. That is, an understanding of fluid, electromagnetic, thermal, and mechanical forces on the micron length scale is necessary in order to understand the operation and function of MEMS devices. In this framework, the methods of actuation and sensing of this new devices have been a critically important topic over the years. There is not a perfect method, and the decision usually depends on the actual device and the specifications of the system. The main actuation and sensing properties used in MEMS are •Piezoresistivity: When a piezoresistive material is stressed, it reduces or increases its ability to transport current. Using this property, movement can be measured as a current difference between the two extremes of a deformed piezoresistive material. 1 Table 1: Comparison between actuation/sensing methods [Kovacs, 1998] Parameter Local circuits DC response Complex Linearity Issues Piezoresistive strain NO YES + +++ High temperature dependance Easy yo integrate Piezoelectric force NO NO ++ ++ High sensitivity Fabrication complex Electrostatic displacement YES YES ++ poor Very simple Low temperature coefficients Thermal strain NO YES + poor Cooling problems Interference with electronics Magnetic displacement NO YES +++ + Very complex Post fabrication Optical displacement NO YES +++ +++ Difficult to implement •Piezoelectricity: Piezoelectric materials deform under the influence of a voltage bias, or reciprocally, under deformation generate a polarization between their extremes. Using this relationship, movement can be controlled or sensed. •Electrostatics: The polarization between two plates generates an electrostatic force between them. This fact can be used to actuate the device. On the other hand, relative movement of two polarized plates generates an induced current that can be sensed, and the movement is proportional to the current. •Thermal: Deformation of the materials due to thermal effects can be used to actuate devices, forcing the increase of temperature in the device. A typical way of achieving the temperature increase is feeding a high current through a conducting material and using the Joule effect. •Electromagnetism: Magnetic fields generated by a current flowing through an spiral can be used to actuate magnetic materials. Similarly, induced current can be used to sense movement of a magnet. •Optics: Reflectivity, transparency, ’admissibility’ of the materials can be used to sense and actuate devices with the help of a light source and a light sensor. The diffraction of the light in a gap, the light patterns of the light through a device, the reflected light in a mirror can be used to extract movement or to induce movement to a MEMS device. Usually, the light source would be a LED or would be carried by an optic fiber. All of them have their advantages and drawbacks (Table 1 and Figure 1), and they are basically related to the selected fabrication method. (See a comparison in [Burns et al., 1995]) Piezoresistive sensing is a common method in engineering to measure strain and displacements. Metal strain gauges are used extensively in engineering. The same principles have been used with semiconductors, and the case apply to doped-silicon or the different layers of material that can be deposited in MEMS (SiO2,Al2O3). Piezoresistive sensing is easy to integrate, and many viable applications exist [Chui et al., 1998], [Tortonese et al., 1993]. However, its temperature dependance and fabrication stresses calibration reduce its market share [Lee, 1997]. Piezoelectric materials are used for actuation and sensing, but the sensing is limited due to their lack of a DC response. Their properties are well known, and have been used for decades. Most of the first sensors used piezoelectric actuation, and it is still used nowadays. However, their high temperature sensitivity, nonlinear working zones and hysteresis prevent from using them more often. When using silicon-based sensors, post-processing is needed to deposit the material. ZnO or PVDF are typical materials used nowadays. Examples of piezoelectric applications could be beam 2 Figure 1: Piezoresistive pressure sensor. Piezoelectric micro-positioner. Analog Devices ADXL150’s electrostatic accelerometer actuation and sensing [Gaucher et al., 1998] and actuation in microscopy [Itoh et al., 1996], [Minne et al., 1995]. Thermal actuation and sensing relays on the use of the thermal deformations of the materials that are used to build the device. The method is easy to implement, and there exist some working devices using this phenomena [Huang and Lee, 2000], [Robert et al., 2003], [Oz and Fedder, 2003]. However, the difficulty of isolating the temperature changes to a fixed area, and the possible interferences with control electronics or other thermally dependent elements, prevents from using this method. [Jonsmann et al., 1999] Magnetic actuation is a common method in the macroworld, however, it is no easily scaled to the MEMS devices. The main problem is the reduction of the achievable forces in a factor of ten thousand when the sizes are reduced by a factor of ten [Niarchos, 2003]. This fact, combined with the constructive difficulties, leaves magnetic actuation application limited. However, successful examples of application exist in the literature, as it could be in gyros [Dauwalter and Ha, 2004], [M Hashimoto and Esashi, 1995] or relays [Tilmans et al., 1999]. Optical actuation and sensing is a desirable method, due to its non-interfering technology. However, although some working devices exist [Lethbridge et al., 1993], [Zook et al., 1995] there are considerable challenges for mass fabrication. The necessity of integrating a light source, building reflecting surfaces and aligning the whole set-up, is time demanding and no batch-fabrication implementation exist. All these problems leave electrostatic actuation and sensing as a really desirable method. Figure 2: Thermal vibromotor [Pai and Tien, 2000]. Optically excited microbeam [Zook et al., 1995]. Micromachined Cu coils [Niarchos, 2003]. 3 Building a capacitor, with the existing fabrication methods is straightforward. One must put together two parallel surfaces and then apply potential difference between the two parts to obtain a good actuator or sensor. This simplicity has made electrostatic actuation and sensing ubiquitous. One can find it in the first MEMS designs to build a gate transistor [Newell, 1968]. Nowadays, capacitive effects are used in resonators [Attia et al., 1998], accelerometers [Kuehnel, 1995], [Brosnihan et al., 1995], optical switches [Juneau et al., 2003], [Sane and Yazdi, 2003], micro-grippers [Chu et al., 1996], micro force gauges [Roessig, 1995], micro-pumps [Teymoori and Abbaspour-Sani, 2002], gyroscopes [Juneau, 1997], [Kranz et al., 2003], pressure sensors [Gupta and Senturia, 1997], RF switches [Huang et al., 2003], and microscopy [Blanc et al., 1996], [Shiba et al., 1998]. Even though practically and economically attractive, capacitive actuation has its own trade-offs and challenges. On-chip amplification is usually needed for capacitive sensing, due to the femptofarad measure that must be achieved. Parasitic capacitances can affect the final read-out. And finally, although large forces can be generated, they can be heavily non-linear. Consequently, the good understanding of the phenomenons that take place is essential to obtain a high performance device with electrostatic actuation and sensing. And this is more relevant given the increasing number of new devices that are continuously designed using these methods of actuation and sensing. 2 Problem Description A basic building block of any electrostatically driven or sensed device is a microbeam. It forms one side of a variable capacity air gap capacitor. Opposite to the microbeam lays the driving or sensing electrode that completes the capacitor. If a voltage is applied to the electrode, a force is generated on the beam that deflects under this action. Alternatively, if the capacitance changes due to a deflection of the microbeam, the charge redistribution and resulting flow of current can be detected. Examples of the typical configurations are shown in Figure 3. (a) (b) (d) (c) V V V V Figure 3: Basic MEMS capacitor configurations (a) Free-end beam. The beam bends under the action of the force. The gap, and consequently the force, is not uniform. Maximum bending at the end. (b) Clamped-clamped beam. The beam bends forming a not uniform gap. Force variable depending on position. Maximum bending in the center. (c) Clamped-clamped beam. A parallel plate added to maintain the gap uniform. Maximum bending in the center which defines the capacitor gap. (d) Guided-end beam. Gap and force uniform. Maximum bending at the extreme of both suspension-beams. 4 When the goal is sensing displacement, a DC polarization voltage is applied to the capacitor, and the generated current is usually detected with a transresistance amplifier [Roessig, 1998]. More sophisticated sensing schemes can also be used to improve the detection. This includes complex electronics designs based on impedance, capacitive or source/drain pick-off [Burstein, 1995]. Charge detection schemes has also been investigated [Seeger and Boser, 2003]. When the goal is driving the beam, an electric load is applied to the microbeam. Depending on the nature of the device, the electric load is composed of a DC polarization voltage and, sometimes, an AC component designed to excite harmonic motions. DC polarization is used to achieve permanent displacements of the beam. Moving optical switches, adjusting elements, closing gate transistors, moving valves or acting micro-grippers are typical applications. However, in most cases, resonant devices are used. In that case, an AC component is added to the driving voltage to excite the harmonic motions of the beam. g wK B M+ _ V Figure 4: Scheme of a parallel plate actuator Figure 4 shows the simplified lumped mass-spring system model of a MEMS device with a parallel plate actuator. To understand the phenomena, one can turn to the energy of the electromechanical system T=1 2M˙ ˆw2;Uk=1 2K ˆw2;Ue=−1 2 ε0Ac (g −ˆw)V2(1) E=T+Uk+Ue(2) where ˆwis the displacement of the moving plate from its initial equilibrium, Tis the kinetic energy of the plate, Ukis the potential energy stored in the spring, Ueis the potential energy stored in the parallel plate capacitor, and Ethe energy of the whole system. The dynamics of the system is derived as follows, using Lagrange’s formulation, d dt µ∂L ∂˙ ˆw¶−∂L ∂ˆw=∂W ∂ˆw(3) being L=T−Uk−Uethe Lagrangian of the system, and introducing the damping force, Fd=−B˙ ˆw as the only contributing force to the work (W) of the system M¨ ˆw+ K ˆw−1 2 ε0Ac (g −ˆw)2V2=−B˙ ˆw(4) 5 This equation is the usual mass-spring-damper equation of dynamics. From this formulation, the force generated between the parallel plates, using basic electrostatics, takes the following form F=1 2 ε0Ac (g −ˆw)2V2(5) where ε0is the dielectric constant, g is the initial gap between the plates, Acis the area of the plates and Vis the applied voltage between the electrodes. As it can be observed, this force is inversely proportional to the gap between the plates of the actuator. As the gap decreases, the generated attractive force increases quadratically. The only opposing force to the electrostatic loading is the mechanical restoring force (K). If the voltage is increased, the gap decreases generating an incremented force. At some point the mechanical forces defined by the spring cannot balance this force anymore. Once reached this state, the electrodes will snap one against the other, and in most cases, the system would be permanently disabled. Consequently, the electrostatic loading has an upper limit beyond which the mechanical force can no longer resist the opposing electrostatic force, thereby leading to the collapse of the structure. This actuation instability phenomenon is known as pull-in, and the associated critical voltage is called the Pull-in Voltage. Several studies have investigated this behavior of microbeams under various loading conditions. The earliest such study may be found in the pioneering work of Nathanson et al. [Nathanson et al., 1967] [Newell, 1968]. In their study of a resonant gate transistor they constructed and analyzed the mass-spring model of electrostatic actuation. They predicted and offered the first theoretical explanation of the so-called pull-in instability. Since then, numerous investigators have analyzed mathematical models of electrostatic actuation in attempts to further understand and control the pull-in instability. Despite more than three decades of work in the area of electrostatically actuated MEMS, the complete dynamics of the electrostatic-elastic system is relatively unexplored. There are a lot of aspects to be clarified. Some studies just center their goal in the immediate application of the sensor, and a simple mass-spring model can approximate the basic dynamics. However, these kind of models cannot predict the inherent nonlinearities of the electrostatic force and the beam deformation ( [Chu et al., 1996], [Casta˜ner and Senturia, 1999]). Other approaches rely on the partial differential equations linearized around the working point. Using this formulation, better results are achieved, but the dynamics only apply for small deflections [Ijntema and Tilmans, 1992]. Other studies analyze the response of a microbeam to a generalized transverse excitation and with axial force using Rayleigh’s energy method to approximate the fundamental natural frequency of the straight, undeflected beam [Tilmans and Legtenberg, 1994]. Recently, some authors have used the nonlinear equation representing the idealized electrostatic structure to analyze the behavior ( [Flores et al., 2003], [Abdel-Rahman et al., 2002]). However, no unified formulation of the problem has been offered. Questions about where, when, and how touchdown occurs are still to be answered. And this knowledge is essential to design and implement the correct control of the new generation of high performance and self-calibrated MEMS devices. 6 3 Model Formulation In this section, the complete idealized model of an electrostatically actuated beam is presented. This model englobes the main characteristics that can be found in a large number of MEMS devices which rely on electrostatic actuation. The analysis of the different participating terms is presented separately, to address each aspect of the dynamics. Finally, the complete formulation is presented together. V L z y x w h b g Figure 5: Basic scheme of a deflected beam 3.1 Mechanical model In MEMS devices, we have a basic structure: the beam. This mechanical component, and its extension, the plate, generate the majority of MEMS sensors and actuators. Consequently, the first step to analyze the behavior of any device is to understand and model the dynamic characteristics of a beam. The deformation of a beam (Figure 5), using the Euler-Bernouilli theory of thin beams [Rao, 1990] is composed of two basic terms [Younis and Nayfeh, 2003], the potential energy generated due to the deformation of the beam Udef =EI 2ZL 0µ∂2ˆw ∂ˆx2¶2 dˆx(6) that it’s proportional to its curvature, ∂2ˆw ∂ˆx2, and the kinetic energy due to its movement T=ρbh 2ZL 0µ∂ˆw ∂ˆ t¶2 dˆx(7) where ˆwis the oscillation amplitude, ρis the density of the beam, b and h are the width and height of the section of the beam, L is the longitude of the beam, E is the Young Modulus and I is the moment of inertia of the cross-section. Typically in MEMS, a beam can also be externally stretched by an axial force ˆ N(ˆ t) (Figure 6). This force could be generated by different sources: thermal load, fabrication stresses, external beam tuning, etc. In this case, another energy term appears that englobes the deformation generated by the external force UN=ˆ N(ˆ t) 2ZL 0µ∂ˆw ∂ˆx¶2 dˆx(8) As can be observed, the deformation is proportional to the axial force. Finally, in the case of large oscillations, the beam movement generates self-stretching forces that actuate as structural damping. This effect can be accounted assuming that an internal force, P, 7 InputForce Figure 6: Vibrating beam oscillating under the influence of an axial force is producing an elongation of the beam. This force would have the following form [Rao and Raju, 2003], [Roessig, 1998] P=bhE 4L ZL 0µ∂ˆw ∂ˆx¶2 dˆx(9) and, substituting this force in (8), we obtain the energy of deformation due to self-stretching [Younis and Nayfeh, 2003] Uint =bhE 8L "ZL 0µ∂ˆw ∂ˆx¶2 dˆx#2 (10) The dynamic equation of the free deflection of an homogeneous beam undergoing bending can be obtained using the Lagrange equations, from the Lagrangian L=T−Udef −UN−Uint (11) and it is written as follows E’I∂4ˆw ∂ˆx4+ρA∂2ˆw ∂ˆ t2−"ˆ N(ˆ t) + E’A 2L ZL 0µ∂ˆw ∂ˆx¶2 dˆx#∂2ˆw ∂ˆx2= 0 (12) where A = bh is the area of the section of the beam, and in this case, the extended Young Modulus, E’ = E/(1 −ν2), is introduced to account for a wide microbeam (plate) where νis the Poisson ratio. For a narrow beam E’ = E. As can be observed, the microbeam dynamics is composed of four terms: the beam resistance to bending, the inertia due to movement, the beam stiffness due to the externally applied axial load and mid-plane stretching due to elongation of the beam. The first three components are treated as linear terms in the equation of motion, whereas the third component is represented by a nonlinear term in the equation of motion. For convenience, and uniformity with other formulations, we introduce the following nondimensional variables w=ˆw g, x =ˆx L, t =ˆ t T(13) where T is a time-scale defined as T = (ρbhL4/(E’I))1/2. Writing down the equation in the nondimensional variables ∂4w ∂x4+∂2w ∂t2−[α1Γ(w, w) + N]∂2w ∂x2= 0 (14) The parameters appearing in equation (14) can be defined as follows α1= 6 ³g h´2, N =ˆ NL2 E’I (15) and the operator Γ is defined as Γ(f1(x, t), f2(x, t)) = Z1 0 ∂f1 ∂x ∂f2 ∂x dx 8 where Pais the static pressure force. As can be observed, the squeeze forces calculation is coupled to the mechanical deflection of the beam [Nayfeh and Younis, 2004]. To approximate the damping forces, one must linearize equation (55) assuming small amplitude motions. This way the gap distance and the pressure of the gap can be expressed as follows d(x, y, t) = g −w(x, y, t) ; P(x, y, t) = Pa+¯ P(x, y, t) (58) where wis the gap reduction and ¯ Pthe pressure variations from the static pressure. Substitution in (55) leads to 12ηeff Pag3µg∂¯ P ∂t −Pa ∂w ∂t ¶=∇2¯ P=∂2¯ P ∂x2+∂2¯ P ∂y2(59) From this equation [Nayfeh and Younis, 2004] has shown that numerical coupled perturbation methods can predict experimental damping forces accurately. If we add the assumption that the capacitor plates are long and narrow (a beam), the equation can be much reduced due to the fact that the fluid movement is only in one direction (y-direction in our device) ∂¯ P ∂t =Pag2 12ηeff ∂2¯ P ∂y2+Pa g ∂w ∂t (60) From this equation, one can solve for ¯ P, obtaining the following force on the capacitors [Senturia, 2001], using Laplace transform Fsq(s) = "96ηeff Lb3 π4g3X nodd 1 n4 1 1 + s αn#sz(s) (61) where αn=g2Pan2π2 12ηeff b2(62) given that z(s) is the input displacement. As we are assuming small amplitudes, the first term of the expansion is a good approximation of the force Fsq(s) = "96ηeff Lb3 π4g3 1 1 + s ωc#sz(s) (63) From this derivation two important parameters arise, the cut-off frequency,ωc ωc=π2g2Pa 12ηeff b2(64) and the squeeze number,σd, σd=π2ω ωc =12ηeff b2 g2Pa ω(65) The squeeze number allow to analyze the behavior of the squeeze film damping forces. When the squeeze number decreases, due to low pressure or low frequencies of oscillation, the fluid force becomes a pure damping force. However, at high frequencies or high squeeze number, a spring force component appears and becomes dominant with the damping force still present. Example of the contributions of each force can be found in [Senturia, 2001]. Similar analysis and discussions are shown by [Andrews et al., 1993] and [Veijola et al., 1995] using the force decomposition derived in [Blech, 1983]. 15 Consequently, squeeze film damping force can be reduced to Fsq =csq(w, σd)∂w ∂t (66) with damping and spring effects depending on σd[Wang et al., 2004]. Finally, the fluid damping effects in the model are the combination of squeeze film and couette film damping, giving a final force Fd=Fsq +Fcou =−ηAov gU+ (P−Pa)·Ac(67) that can be generalized as Fd= (csq +ccou)∂w ∂t = ˆcd ∂w ∂t (68) 3.4 Lumped system The complete set of equations defining the behavior of the system can be obtained linking the different energies and non-conservative forces acting in the system. The kinetic energy is defined in (7) T=ρbh 2ZL 0µ∂ˆw ∂ˆ t¶2 dˆx(69) The potential energy is composed of mechanical (6),(8),(10) and electrostatic terms (30) U=EI 2ZL 0µ∂2ˆw ∂ˆx2¶2 dˆx+ˆ N(ˆ t) 2ZL 0µ∂ˆw ∂ˆx¶2 dˆx+bhE 8L "ZL 0µ∂ˆw ∂ˆx¶2 dˆx#2 +εV 2 2Zv|∇ψ|2dv(70) The fluid damping is the only non-conservative force (67) Fd=−ηAov gU+ (P−Pa)·Ac(71) Consequently, using Lagrange formulation and non-dimensional variables, the dynamics of the system is as follows: ∂4w ∂x4+∂2w ∂t2−[α1Γ(w, w) + N]∂2w ∂x2=γV 2|∇ψ|2−12L4 E’ h3T·−ηAov gU+ (P−Pa)·Ac¸(72) given that the electrostatic potential and the fluid pressure satisfy the following conditions ²2µ∂2ψ ∂x2+ a2∂2ψ ∂y2¶+∂2ψ ∂z2= 0 (73) 12ηeff ∂Pd ∂t =∇[d3P∇P](74) Linking the different formulations previously derived , the dynamics of the system can be reduced to [Abdel-Rahman et al., 2003]: ∂2w ∂t2+c∂w ∂t +∂4w ∂x4−[α1Γ(w, w) + N]∂2w ∂x2=γV 2|∇ψ|2(75) 16 w(0, t) = w(1, t) = 0, w0(0, t) = w0(1, t) = 0 And the parameters appearing in equation (75) can be defined as follows c=ˆcdL4 E’ I T , N =ˆ NL2 E’I α1= 6 ³g h´2, γ =6ε0L4 E’ h3g(76) Equation (75) translates to the following formulation once the electrostatic force is approximated ∂2w ∂t2+c∂w ∂t +∂4w ∂x4−[α1Γ(w, w) + N]∂2w ∂x2=κV2 (1 −w)2(77) where κ=6Cnε0L4 E’ h3g3using fringing fields correction. 4 Model solution Once the model has been derived, one can analyze the behavior of the system. In this section the different approaches to understand the system are presented and formulated. With each approach the advantages and problems are presented, as well as, the implications to the stability of the system. 4.1 Static solution In the case of searching for the static solution of the system, the time-derivatives of the system must be set to zero. Under these premises, only potential energy terms remain in our system and the static solutions correspond to the equilibrium positions of the potential energy of the system, that is dU dˆw= 0 (78) Consequently, the static deformation wsof the beam, under the action of a electrostatic forcing Vp can be calculated from equation (75), if the time-derivatives are set to zero. This way, the remaining terms are only position-dependant, and the partial differential equations disappear: d4ws dx4−[α1Γ(ws, ws) + N]d2ws dx2=γV 2 p|∇ψ|2(79) ws= 0 and dws dx= 0 at x = 0 and x = 1 (80) Unfortunately, equation (79) do not generate a closed-form solution, due to its implicit nature. For this reason, numerical methods must be used to solve the problem. A possibility is shooting methods combined with nonlinear boundary-value problem solution as in [Abdel-Rahman et al., 2003]. They apply the method to this model without fringing-fields correction (Cn= 1). d4ws dx4−[α1Γ(ws, ws) + N]d2ws dx2=κV2 (1 −w)2(81) They show good agreement to experimental results, and argue that the inclusion of internal stretching is essential to predict real displacements. This approach allows to numerically calculate 17 the exact Static Pull-in Voltage using the same numerical method. Their analysis shows that neglecting the nonlinear effects leads to underestimating the stability limits of the system. The travel range taking into account the nonlinearities can be doubled. Another option is to ignore the internal stretching, what reduces the model complexity, and use a method as the backward Euler algorithm to solve for the static displacement as in [Ijntema and Tilmans, 1992]: d4ws dx4−Nd2ws dx2=κV2 (1 −w)2(82) [Tilmans and Legtenberg, 1994] solved the same static problem using the Rayleigh-Ritz method assuming a combination of trial functions. They used this formulation to generate an analytical expression for the pull-in voltage, based on energy methods. Even with the needed approximations to solve the equations, the calculated values of the pull-in voltage were in good agreement with the results of experiments they conducted on resonators of various lengths. The system approximation generates good results while large amplitudes are not taken into account. In [Zhou and Yang, 2003] numerical solutions are shown using the equation (82) and finite elements analysis. General numerical solutions using finite elements with reduced-order energy equations are presented in [Elata et al., 2003] using relaxation techniques. FEM solutions allow to handle the complete deformation of the device, without focusing on the maximum amplitude, but large computational time is needed. 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 16 14 12 10 8 6 4 2 0 2 4 x 10 12 normalized displacement (y/g0) Energy Evolution energy profile with increasing voltage 10 V 30 V 40 V 55.43V 60.34V 65V ypin=g0/3 Stable and unstable equilibrium points for VDC=30V Initial energy and energy at unstable equilibrium point for DPV=55.43V Figure 9: The Potential Energy levels of the a parallel-plate actuator system depend on the position relative to the gap. Energy of the system versus normalized displacement for different applied voltages are displayed for an example, including the Static Pull-in Voltage (60.34 V) and the Dynamic Pull-in Voltage (55.43 V). The stable equilibrium corresponds with the static displacement of the device. Some elaborated solutions and behavior analysis are derived in [Bernstein et al., 2000] and [Pelesko, 2001a] directly form the differential equations. To arrive to the solutions, a simplified membrane model is used where the plate inertial and bending effects are neglected. However, numerical implicit formula solution is also needed to evaluate the static solution. d2u dx2=β (1 + u)2(83) 18 This analysis allows to define stability conditions based on implicit eigenvalue equations. Most authors work with the mass-spring-damper model, as in Figure 4. This model losses insight on the complete behavior of the system, but allow to analyze the system analytically, producing important information for the design process. The behavior of the beam can be approximated to that of a non-linear spring for a given deformation mode, as it has been shown in (19), and approximations can also be obtained for the electrostatic force and damping, giving way to the following formulation Meff,i ·¨ ˆqi+Ceff,i ·˙ ˆqi+Keff,i ·ˆqi+K3,eff,i ·ˆq3 i=Fe(84) In the static case, (84) simplifies to Keff,i ·ˆqi+K3,eff,i ·ˆq3 i=Fe(85) being a non-linear mass-spring equation. This model characterizes the beam stiffening due to large deformations, that reduces the effective travel range of the beam [Roessig, 1998]. However, in most cases small amplitude of oscillation is considered [Vinokur, 2002], [Gretillat et al., 1997] , allowing to use the linear formulation K·ˆqi=−1 2 ε0AcV2 g2(1 −w)2(86) With this model, the classical Static Pull-in Voltage (SPV ) equation is obtained (Figure 9), SPV =s8 27 K g3 0 ε0A;ypin =g0 3(87) which indicates the maximum voltage that can be applied without getting snapping. Substitution of the voltage in the dynamics equation gives the maximum travel range in the static case, which is one-third of the initial gap [Senturia, 2001]. 4.2 Dynamic solution If we want to analyze the transient response of the microbeam when a variable voltage load V(t) is applied, the complete evolution of the energy of the system has to be taken into account. To obtain solutions, the full set of equations (75) must be used. Typical cases where the transient is of interest include micro-switches or mirror positioning, where the time response is of great interest. In this cases, the voltage is usually applied as a step-function or a ramp-function. The use of numerical simulation to obtain the behavior of the system is mandatory if the complete set of equations is used. MEMS exhibit non-linearities even when the displacements are small, and the complete equations are needed to capture all the behavioral aspects. The main dynamic nonlinear effects that will be detected on parallel-plate actuated MEMS are the following [Rand, 2003]: •Spring stiffening: The effect appears due to large amplitudes of oscillation. The deformation of the beam cannot be considered linear anymore, and increases the beam resistance to deformation (19). The resulting non-linear equation corresponds to the Duffing equation (Figure 10a). 19 •Spring softening: The electrostatic force function (36) can be approximated using Taylor series. In that case, if only the first and second term are used, a negative spring term appears in the system equations. This fact is usually detected as a natural frequency reduction while increasing the voltage bias. This fact is used in some cases to adjust and trim the frequency of MEMS resonators [Painter and Shkel, 2003]. •Parametric excitation: The spring softening generated by the electrostatic force can derive to parametric excitation when an oscillatory force is used. In that case, the system behavior is governed by a Mathieu equation [Butikov, 2004]. Particular analysis can be carried out to analyze the parametric resonances and instabilities (Figure 10b). •Hysteresis: Associated to the Duffing nonlinearity (Figure 10a), the system can derive to have bifurcation points that generate hysteresis regions in the behavior of the system [Gui et al., 1998], [Kaajakari et al., 2005]. •Chaos regions: Some works have analyzed the nonlinear behavior of parallel-plate actuated MEMS detecting existence of chaotic regions that could restrict the stable range of actuation of the devices [Liu et al., 2004] [Bienstman et al., 1998] [Wang et al., 1998]. (b) (a) Figure 10: (a) Characteristic non-linear Duffing equation behavior of the frequency response of a parallel plate oscillator. As the amplitude increases, a hysteresis appears [Gui et al., 1998] ; (b) Characteristic profile of exponential growth during parametric excitation [Napoli et al., 2004] A complete simulation of the system is presented in [Nayfeh and Younis, 2004], where the modeling and simulation under the effect of squeeze-film damping is analyzed. They use the compressible Reynolds equation coupled with the equation governing the plate deflection (72-74). The model accounts for the electrostatic forcing of the capacitor air-gap, the restoring force of the microplate and the applied in-plane loads. Perturbation methods are used to derive an analytical expression for the pressure distribution. This expression is then substituted into the plate equation, which is solved in turn using a finite-element method for the structural mode shapes, the pressure distributions, the natural frequencies and the quality factors. Without taking the damping into account, some works analyze the electro-mechanical behavior. Analysis of the equations is carried out in [Xie et al., 2003] using a nonlinear modal analysis approach based on the invariant manifold method. Using Galerkin method, the nonlinear partial differential governing equation is decoupled into a set of nonlinear ordinary differential equations. 20 Then the invariant manifold method is used to obtain the associated nonlinear modal shapes, and modal motion governing equations. The model allows to examine the nonlinearities and the pullin phenomena. Similar results using shooting methods combined with nonlinear boundary-value problem where presented in [Abdel-Rahman et al., 2002]. Using a simplified plate model (83), in [Flores et al., 2003] they obtain solutions of the system operated in viscous regime. This simplified mathematical model allows to study a parabolic equation of reaction-diffusion type. A central result of the paper is that when the applied voltage is beyond the critical voltage where steady-state solutions cease to exist, the solution touches down in finite time. Bounds on the touchdown time are computed and the structure of solutions near touchdown are investigated. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 4.62 4.6 4.58 4.56 4.54 4.52 4.5 4.48 4.46 4.44 x 10 13 normalized displacement (y/g0) Energy Energy evolution Potential energy bound System energy evolution Initial energy level of the system Stable equilibrium position 10210 1 100101102103104105 55 56 57 58 59 60 61 Quality Factor Pull-in Voltage (V) Dynamic Pull-in Voltage Evolution Static Pull-in Voltage = 60.34 V Dynamic Pull-in Voltage = 55.43 V Q=1.2 Q=1000 (a) (b) Figure 11: (a) Evolution of system’s energy of an example when a 30 V step-function is applied. The Quality Factor in the example is 30. The initial energy corresponds to the potential energy (mechanical and electrostatic). When the motion begins, the potential energy is converted to kinetic energy and dissipation due to damping forces. The system’s energy descends until reaching the stable equilibrium position; (b) Evolution of the pull-in voltage as a function of the Quality Factor in a example. For high-Q environments, the pull-in voltage corresponds to the Dynamic Pull-in Voltage. For low-Q environments, the pull-in voltage corresponds to the Static Pull-in Voltage. The calculation were done with a linear mass-spring-damper model. The equations can be reduced to mass-spring-damper formulation as in [Krylov et al., 2005] and [Krylov and Maimon, 2004] in order to highlight leading dynamical phenomena through analysis of simplified expressions. They develop a model using the Galerkin procedure with normal modes as a basis. It accounts for the distributed nonlinear electrostatic forces, nonlinear squeezed film damping, and rotational inertia of a mass carried by the beam. Special attention is paid to the dynamics of the beam near instability points. The results generated by the model, and confirmed experimentally, show that nonlinear damping leads to shrinkage of the spatial region where stable motion is realizable. The model is useful to generate conclusions about the stability using the simplest model of a parametrically excited system described by Mathieu and Hills equations. Energy methods are used in [Ligterink et al., 2005] to analyze the transient behavior between pull-in and release states. The concept of dynamic pull-in is addressed as well as hysteresis phenomena. No evolution analysis are performed. Using mass-spring-damper models, linked to FEM analysis, interesting results can also be obtained. In [Han et al., 2005], model order reduction techniques are used to reduce the transient 21 analysis time. To do this, an open-source software performs model order reductions via the block Arnoldi algorithm directly to ANSYS finite element models. On the other hand, direct analysis over the linear model can be useful in several applications (Figure 11a). The nonlinear behavior of the system with a simple mass-spring-damper model is analyzed in [ZHAO et al., 2005], [Casta˜ner et al., 1999], [Minami et al., 1999]. In [Gupta and Senturia, 1997], they used the system analysis to predict pull-in times and derive the Dynamic Pull-in Voltage (DPV) yuns =g0 2;DPV =s1 4 K g3 0 ε0A(88) which indicates the maximum voltage that can be applied as a step-function to the system without producing snapping in vacuum environment. A extended discussion on energy-dependence of the Dynamic Pull-in Voltage (Figure 11b) can be found in [Varghese et al., 1997] and [Fargas-Marques, 2001]. 4.3 Oscillatory solution In multiple applications in MEMS sensors and actuators, the system is oscillated at a fixed frequency. An alternating voltage is applied to the system to maintain the oscillation. The case of oscillatory load is, then, a sub-case of the dynamic solution. To analyze the oscillatory case, the transient response is neglected and the efforts are concentrated on the stationary oscillation. The microbeam deformation under an electrostatic excitation (V(t) = V p +v(t)) is composed of a static component (ws(x)) and a dynamic component (u(x, t)), due to the AC forcing voltage: w(x, t) = ws(x) + u(x, t) (89) To solve the oscillatory case, we substitute (89) in the dynamic equation of the system (75) and obtain ∂2(ws(x)+u(x,t)) ∂t2+c∂(ws(x)+u(x,t)) ∂t +∂4(ws(x)+u(x,t)) ∂x4 −[α1Γ((ws(x) + u(x, t)),(ws(x) + u(x, t))) + N]∂2(ws(x)+u(x,t)) ∂x2 =γ(Vp+v(t))2|∇ψ|2(90) The equation can be simplified eliminating the expressions with null terms, and it turns to ∂2u ∂t2+c∂u ∂t +∂4ws ∂x4+∂4u ∂x4 −[α1(Γ(ws, ws) + 2Γ(ws, u) + Γ(u, u)) + N] (∂2ws ∂x2+∂2u ∂x2) =γ(Vp+v(t))2|∇ψ|2(91) Once in this point, to develop the equation further, a possibility is to approximate the electrostatic force, assuming no fringing fields, using equation (32). This way the oscillation is defined by ∂2u ∂t2+c∂u ∂t +∂4ws ∂x4+∂4u ∂x4 −[α1(Γ(ws, ws) + 2Γ(ws, u) + Γ(u, u)) + N] (∂2ws ∂x2+∂2u ∂x2) =α2(Vp+v(t))2 (1−(ws+u))2(92) where now α2=γ/g2=6ε0L4 E’ h3g3. 22 The formulation can be much reduced if the electrostatic force is expanded in Taylor series around the equilibrium position ∂2u ∂t2+c∂u ∂t +∂4ws ∂x4+∂4u ∂x4 −[α1(Γ(ws, ws) + 2Γ(ws, u) + Γ(u, u)) + N] (∂2ws ∂x2+∂2u ∂x2) =α2(V2 p+ 2Vpv(t) + v(t)2)³1 β2+2 β3u+3 β4u2+4 β5u3+O(u4)´(93) and then, rearranging terms, the static solution (79) can be eliminated and only the oscillating solution is conserved, simplifying the equation to ∂2u ∂t2+c∂u ∂t +∂4u ∂x4 −α1[2Γ(ws, u) + Γ(u, u)] ∂2ws ∂x2 −[α1(Γ(ws, ws) + 2Γ(ws, u) + Γ(u, u)) + N]∂2u ∂x2 =α2V2 p³2 β3u+3 β4u2+4 β5u3+O(u4)´ +2α2Vpv(t)³1 β2+2 β3u+3 β4u2+4 β5u3+O(u4)´ +α2v(t)2³1 β2+2 β3u+3 β4u2+4 β5u3+O(u4)´(94) Depending on the driving voltages and the accepted error, the final equation can be selected. However, numerical simulation will be needed to obtain evolution results. Complete simulations based on the theoretical framework exist in the literature. In [Abdel- Rahman et al., 2003], shooting methods combined with nonlinear boundary-value problem are used to solve the existing eigenvalue problem. The vibrations around the deflected position of the microbeam are solved numerically for various parameters to obtain the natural frequencies and mode shapes. The results are compared with experimental results available in the literature with good agreement. In [Nayfeh and Younis, 2004], perturbation methods are used to obtain the mode shapes and frequencies including the coupled effects of the squeeze-film damping that were just approximated in the previous analysis. Following the same study, in [Younis et al., 2004] they present a methodology to simulate the transient and steady-state dynamics of microbeams undergoing small or large motions actuated by combined DC and AC loads. They use the model to produce results showing the effect of varying the DC bias, the damping, and the AC excitation amplitude on the frequency-response curves. In their analysis they detect the existence of dynamic effects that can produce pull-in with electric loads much lower than that predicted based on static analysis. In [Ijntema and Tilmans, 1992], the dynamic behavior is modeled using energy methods to obtain a spring-mass-damper model. The fundamental frequency is approximated using Rayleigh’s energy method where the microbeam motion is linearized around the deflected shape obtained as a solution of the static problem. The work is extended in [Tilmans and Legtenberg, 1994] where the fundamental natural frequency obtained from Rayleigh’s energy method is compared to the experimentally obtained fundamental natural frequency. They found out that the results obtained from the expression were only valid for small dc polarization voltages away from the pull-in voltage. Their method takes into account the axial load and large amplitude effects. The mass-spring-damper formulation obtained via Galerkin procedure in [Krylov et al., 2005] allows to study the parametric resonance behavior of the system. They show that parametric stabilization can be obtained. The model summarizes the main nonlinearities for a given frequency. Similar analysis are carried out with a parametric model in [Napoli et al., 2004] They show that the underlying linearized dynamics of the system are those of a periodic system described by a Mathieu 23 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 12 Normalized displacement (y/g0) Energy Energy evolution Maximum stable amplitude of oscillation Achieved amplitude of oscillation Potential energy upper bound VDC-VAC Potential energy lower bound VDC+VAC Initial position Unstable equilibrium of VDC+VAC curve Total energy evolution Stable permanent oscillation loop (a) Initial energy 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 12 Normalized displacement (y/g0) Energy Energy evolution Maximum stable amplitude of oscillation Potential energy upper bound VDC-VAC Potential energy lower bound VDC+VAC Initial position Unstable equilibrium of VDC+VAC curve Total energy evolution Unstable oscillation loop (b) Snapping of electrodes Initial energy Figure 12: The potential energy curves bound the system oscillation. In (a) a stable oscillation is obtained for the example system with 19 VDC bias voltage and a 7 VAC amplitude while in (b) the oscillation is unstable with 20 VDC and 7 VAC. Beginning from the static initial position, the amplitude of oscillation increases until it reaches the unstable equilibrium point at VDC +VAC, resulting in snapping. equation. Experimental results confirm the validity of the model, and in particular, illustrate that parametric resonance phenomena occur in capacitively actuated micro-cantilevers. Finite element approaches are also valid to obtain the behavior of the system [Gretillat et al., 1997]. However computation times can be quite demanding in the case of non-linear coupling. In [Hung, 1997], it is shown that a way of solving the simulation of the system is rewriting the solution as a sum of orthogonal basis functions, that correspond to the oscillation modes. They show the feasibility using an initial model with internal tension and damping. The obtained low-order models are quicker for numerical modeling. Multiple authors have modeled a microbeam under electrostatic actuation as a single-degree-of- freedom spring-mass-damper system. The model assumes a linear spring, thus neglects midplane stretching effects. They use this model to generate an analytical expression for the fundamental natural frequency as a function of the dc polarization voltage. Both this expression and the experiments they carry on a resonator show that increasing the dc polarization voltage decreases the fundamental natural frequency [Vinokur, 2002], [Seeger, 1997], [Sung et al., 2003]. Using a spring-mass-damper model and energy methods, the AC Pull-in Voltage is presented and analyzed in [Fargas-Marques, 2001] as the combination of VDC and VAC that can lead the system to snapping (Figure 12). Numerical and experimental results are presented to validate the concept. Similar results are presented in [Seeger and Boser, 2002] for double-sided actuated oscillators near the mechanical resonant frequency and amplitudes comparable to the actuator gap. They show that at resonance, the structure can move beyond the well-known pull-in-limit but is instead limited to 56% of the gap by resonant pull-in. Above the resonant frequency, the structure is not limited by pull-in and can theoretically oscillate across the entire gap. 24 [Niarchos, 2003] Niarchos, D. (2003). Magnetic mems: key issues and some applications. Sensors and Actuators A. [Nishiyama and Nakamura, 1990] Nishiyama, H. and Nakamura, M. (1990). Capacitance of a strip capacitor. Components, Hybrids, and Manufacturing Technology, IEEE Transactions on. [Oz and Fedder, 2003] Oz, A. and Fedder, G. (2003). Rf cmos-mems capacitor having large tuning range. TRANSDUCERS, Solid-State Sensors, Actuators and Microsystems, 12th International Conference on,. [Pai and Tien, 2000] Pai, M. and Tien, N. C. (2000). Low voltage electrothermal vibromotor for silicon optical bench applications. Sensors and Actuators. [Painter and Shkel, 2003] Painter, C. C. and Shkel, A. M. (2003). Active structural error suppression in mems vibratory rate integrating gyroscopes. IEEE SENSORS JOURNAL. [Pelesko, 2001a] Pelesko, J. (2001a). Electrostatic field aproximations and implications for mems devices. Proceedings of ESA. [Pelesko, 2001b] Pelesko, J. (2001b). Multiple solutions in electrostatic mems. Proceedings of Modeling and Simulation of Microsystems 2001, Hilton Head. [Pelesko and Triolo, 2000] Pelesko, J. and Triolo, A. (2000). Nonlocal problems in mems device control. In Technical Proceedings of the 2000 International Conference on Modeling and Simulation of Microsystems. [Pelesko and Triolo, 2001] Pelesko, J. and Triolo, A. (2001). Nonlocal problems in mems device control. Journal of Engineering Mathematics, 41. [Pelesko and Bernstein, 2003] Pelesko, J. A. and Bernstein, D. H. (2003). Modeling MEMS and NEMS. Chapman & Hall/CRC. ISBN: 1-58488-306-5. [Rand, 2003] Rand, R. H. (2003). Lecture Notes on Nonlinear Vibrations, version 45. Available online at http://www.tam.cornell.edu/randdocs/. [Rao and Raju, 2003] Rao, G. and Raju, K. (2003). Large amplitude free vibrations of beams - an energy approach. ZAMM - Journal of Applied Mathematics and Mechanics, 83(7):493 – 498. [Rao, 1990] Rao, S. (1990). Mechanical Vibrations. Addison-Wesley, 2nd edition. [Robert et al., 2003] Robert, P., Saias, D., Billard, C., Boret, S., Sillon, N., Maeder-Pachurka, C., Charvet, P., Bouche, G., Ancey, P., and Berruyer, P. (2003). Integrated rf-mems switch based on a combination of thermal and electrostatic actuation. TRANSDUCERS, Solid-State Sensors, Actuators and Microsystems, 12th International Conference on,. [Rocha et al., 2004] Rocha, L., Cretu, E., and Wolffenbuttel, R. (2004). Pull-in dynamics: analysis and modeling of the transitional regime. In Micro Electro Mechanical Systems,. [Roessig, 1995] Roessig, T. A. W. (1995). Surface micromachined resonant force transducers. Master’s thesis, U.C. Berkeley. [Roessig, 1998] Roessig, T. A. W. (1998). Integrated MEMS Tuning Fork Oscillators for Sensor Applications. PhD thesis, U.C. Berkeley. 31 [Sane and Yazdi, 2003] Sane, H. and Yazdi, N.and Mastrangelo, C. (2003). Application of sliding mode control to electrostatically actuated two-axis gimbaled micromirrors. In American Control Conference, 2003. Proceedings of the 2003. [Saucedo-Flores et al., 2004] Saucedo-Flores, E., Ruelas, R., Flores, M., Ying, C., and Jung-chih., C. (2004). Dynamic behavior modeling of mems parallel plate capacitors. PLANS 2004. Position Location and Navigation Symposium (IEEE Cat. No.04CH37556).IEEE . Piscataway, NJ, USA, pages pp.15–19. [Seeger and Boser, 2003] Seeger, J. and Boser, B. (2003). Charge control of parallel-plate, electrostatic actuators and the tip-in instability. Microelectromechanical Systems, Journal of. [Seeger and Boser, 2002] Seeger, J. I. and Boser, B. E. (2002). Parallel-plate driven oscillations and resonant pull-in. In Solid-State Sensor, Actuator and Microsystems Workshop Hilton Head Island. [Seeger, 1997] Seeger, J.I. Crary, S. (1997). Stabilization of electrostatically actuated mechanical devices. Solid State Sensors and Actuators, 1997. TRANSDUCERS ’97. [Senturia, 2001] Senturia, S. (2001). Microsystem Design. Kluwer Acadenic Publishers, 1st edition. [Shiba et al., 1998] Shiba, Y., Ono, T., Minami, K., and Esashi, M. (1998). Capacitive afm probe for high speed imaging. ns. of the IEE of Japan, 118 - E(12):647–650. [Sung et al., 2000] Sung, S., Lee, J., Kang, T., and Song, J. W. (2000). Development of a tunable resonant accelerometer with self-sustained oscillation loop. In IEEE National Aerospace and Electronics Conference. [Sung et al., 2003] Sung, S., Lee, J. G., and Kang, T. (2003). Development and test of mems accelerometer with self-sustatined oscillation loop. Sensors and Actuators A. [Teymoori and Abbaspour-Sani, 2002] Teymoori, M. and Abbaspour-Sani, E. (2002). A novel electrostatic micromachined pump for drug delivery systems. In Semiconductor Electronics, 2002. Proceedings. ICSE 2002. IEEE International Conference on. [Tilmans et al., 1999] Tilmans, H., Fullin, E., Ziad, H., van der Peer, M., Kesters, J., and al. (1999). A fully-packaged electromagnetic microrelay. Micro Electro Mechanical Systems. [Tilmans and Legtenberg, 1994] Tilmans, H. A. C. and Legtenberg, R. (1994). Electrostatically driven vacuum-encapsulated polysilicon resonators part ii. theory and performance. Sensors and Actuators A: Physical. [Tortonese et al., 1993] Tortonese, M., Barrett, R., and Quate, C. (1993). Atomic resolution with an atomic force microscope using piezoresistive detection. Appl. Phys. Lett, 62:8340–8363. [Varghese et al., 1997] Varghese, M., Amantea, R., Sauer, D., and Senturia, S. D. (1997). Resistive damping of pulse-sensed capacitive position sensors. In TRANSDUCERS ’97, International Conference on Solid-state Sensors and Actuators. [Veijola et al., 1995] Veijola, T., Kuisma, H., Lahdenper, J., and Ryhnen, T. (1995). Equivalentcircuit model of the squeezed gas film in a silicon accelerometer. Sensors and Actuators A: Physical. 32 [Veijola and Turowski, 2001] Veijola, T. and Turowski, M. (2001). Compact damping models for laterally moving microstructures with gas-rarefaction effects. Microelectromechanical Systems. [Vinokur, 2002] Vinokur, R. Y. (2002). Feasible analytical solutions for electrostatic parallel-plate actuator or sensor. Journal of Vibration and Control. [Wang et al., 2004] Wang, X., Liu, Y., Wang, M., and Chen, X. (2004). The effect of air damping on the planar mems structures. In HDP’04. [Wang et al., 1998] Wang, Y. C., Adams, S. G., Thorp, J. S., MacDonald, N. C., Hartwell, P., and Bertsch, F. (1998). Chaos in mems, parameter estimation and its potential application. IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS. [Xie et al., 2003] Xie, W. C., Lee, H. P., and Lim, S. P. (2003). Nonlinear dynamic analysis of mems switches by nonlinear modal analysis. Nonlinear Dynamics, 31:243256. [Yang and Senturia, 1996] Yang, Y. and Senturia, S. (1996). Numerical simulations of compressible squeezed-film damping. In Solid-State Sensors and Actuators Workshop, Late News Session. [Yang et al., 1997] Yang, Y.-J., Gretillat, M.-A., and Senturia, S. (1997). Effect of air damping on the dynamics of nonuniform deformations of microstructures. Solid State Sensors and Actuators, 1997. TRANSDUCERS ’97. [YEH et al., 2001] YEH, B. Y., LIANG, Y. C., and TAY, F. E. H. (2001). Mathematical modelling on the quadrature error of low-rate microgyroscope for aerospace applications. Analog Integrated Circuits and Signal Processing. [Younis et al., 2004] Younis, M. I., Abdel-Rahman, E., and Nayfeh, A. H. (2004). Global dynamics of mems resonators under superharmonic excitation. In Proceedings of the 2004 International Conference on MEMS, NANO and Smart Systems. [Younis and Nayfeh, 2003] Younis, M. I. and Nayfeh, A. H. (2003). A study of the nonlinear response of a resonant microbeam to an electric actuation. Nonlinear Dynamics, 31:91117. [Zhang et al., 2002] Zhang, W., Baskaran, R., and Turner, K. L. (2002). Effect of cubic nonlinearity on auto-parametrically amplified resonant mems mass sensor. Sensors and Actuators A. [ZHAO et al., 2005] ZHAO, X., REDD, C. K., and NAYFEH, A. H. (2005). Nonlinear dynamics of an electrically driven impact microactuator. Nonlinear Dynamics. [Zhou and Yang, 2003] Zhou, Y.-H. and Yang, X. (2003). Numerical analysis on snapping induced by electromechanical interaction of shuffling actuator with nonlinear plate. Computers and Structures. [Zook et al., 1995] Zook, J., Burns, D., Herb, W., Guckel, H., Kang, J., and Ahn, Y. (1995). Optically excited self-resonant strain transducers. In Transducers’95, 8th International Conference on Solid-State Sensors and Actuators, volume 2. 33