# Dynamic Analysis of Multilayers Based MEMS Resonators

Authors: Hassen M. Ouakad; Abdulrahman M. Alofi; Ali H. Nayfeh
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The dynamic behavior of a microelectromechanical system (MEMS) parallel and electrically coupled double-layers (microbeams) based resonator is investigated. Two numerical methods were used to solve the dynamical problem: the reduced-order modeling (ROM) and the perturbation method. The ROM was derived using the so-called Galerkin expansion with considering the linear undamped mode shapes of straight beam as the basis functions. The perturbation method was generated using the method of multiple scales by direct attack of the equations of motion. Dynamic analyses, assuming the above two numerical methods were performed, and a comparison of the results showed good agreement. Finally, a parametric study was performed using the perturbation on different parameters and the results revealed different interesting features, which hopefully can be useful for some MEMS based applications.


---

## Body

## 1. Introduction

Microelectromechanical systems (MEMS) were mostly developed during the industrial revolution in the late of 20th century [1]. From that moment, scientists have continuously investigated the potential of these devices to apply them for different applications. Nowadays, their existence became indispensable in many fields [2, 3] and the continuous development hopefully will lead to more importance in future.When dealing with MEMS based devices, interesting dynamical behaviors can be arising due to the nonlinearity of those devices [4–6]. Among them is their hysteresis (softening and hardening) behavior during the resonance phenomenon. The hardening effect shifts the linear resonance profile to the right and hence produces higher frequencies for the higher amplitude values [7]. A main source of this effect for the clamped-clamped microbeam is the mid-plane stretching. On the other hand, the opposite can be said for the softening effect, which depends on the quadratic nonlinearity effect that may arise from the nonlinear electrical forces [8].Solving nonlinear dynamical behaviors for MEMS devices is fundamental, since it helps in accurately characterizing and designing them to obtain the desired features quickly and effectively. However, the resolution could be sometime cumbersome and many researchers have struggled to find an effective way to tackle this problem [9–12]. The key techniques in solving such nonlinear problems involve reducing the order of the partial differential governing equation, which sometime is difficult to be solved, into ordinary differential equations, which are much easier to deal with. Different approaches are suggested to solve equation; the reduced-order modeling (ROM) and the perturbation analysis are among those approaches.In recent years, many researchers extensively used reduced-order modeling (ROM) techniques to obtain the structural behavior of MEMS microstructures [13, 14]. It became nowadays among the most common numerical methods used in the MEMS community. In 2003, Younis et al. [15] presented a new ROM (macromodel) to solve for the nonlinear partial differential equation. They found that their method is more attractive than the finite element method in terms of accuracy and cost. Nayfeh et al. [16] reviewed the work of ROM in MEMS devices. They classified the ROM into two main categories: node and domain methods.Perturbation method is also another approach used in the investigation. Turner and Andrews [17] used the perturbation method to obtain an approximation for the nonlinear resonance frequency of a microbeam. They modeled the problem by using a spring mass model, and they included a cubic restoring force to represent the mid-plane stretching. Younis and Nayfeh [6] used the method of multiple scales, as a perturbation method, in a model that accounted for the electrical loads (both DC and AC) and for the mid-plane stretching for a clamped-clamped single-microbeam resonator.In a previous study [18], we examined the static behavior of a double-layers based clamped-clamped microbeams MEMS based actuator configuration. This work considers the nonlinear dynamical behaviors for this MEMS structure by using the ROM and the perturbation method. This is due to the different features this device has as it is compared with the single microbeam. For example, the double-microbeam configuration proved its capability in reducing the actuation and the pull-in voltage, and nowadays there are many available investigations that used this type of configuration [19–23].We will start our present numerical investigation by discussing first the problem formulation. Then, the reduced-order model for double-layers MEMS based resonator will be derived and followed by the analysis for the perturbation method. Afterward, a comparison will be carried out for each method to validate both approaches. Finally, a parametric study will be done to investigate some dynamical features for the microactuator.


## 2. Problem Formulation

The proposed model is shown in Figure1. The suggested MEMS device is mainly made of two parallel electrically actuated microbeams. Both microbeams are clamped-clamped and actuated by both a static DC voltage (of amplitude VDC) and a harmonic AC bias (of amplitude VAC).Figure 1
3D schematic of a MEMS double-layers MEMS based resonator configuration.The governing equations of motion for upper and lower microbeams of the MEMS multilayers based resonator are given, respectively, as follows:(1)EI∂4w^1∂x^4+ρbh∂2w^1∂t^2+c^∂w^1∂t^=EA2L∫0L∂w^1∂x^2dx^∂2w^1∂x^2+ε0bVDC+VACcosΩ^t^22d1-w^12-ε0bVDC+VACcosΩ^t^22d2+w^1-w^22,(2)EI∂4w^2∂x^4+ρbh∂2w^2∂t^2+c^∂w^2∂t^=EA2L∫0L∂w^2∂x^2dx^∂2w^2∂x^2+ε0bVDC+VACcos⁡Ω^t^22d2+w^1-w^22and the associated clamped-clamped boundary conditions are(3)w^10,t^=0,w^1L,t^=0,∂w^1∂x^0,t^=0,∂w^1∂x^L,t^=0,w^20,t^=0,w^2L,t^=0,∂w^2∂x^0,t^=0,∂w^2∂x^L,t^=0.The functions w^1x^,t^ and w^2x^,t^ are, respectively, the in-plane deflections of the upper and lower microbeam, respectively. For convenience, the equation of the system will be normalized while considering the following nondimensional parameters:(4)w1=w^1d1,w2=w^2d1,x=x^L,t=t^T,where T is a time scale parameter and is chosen to be ρbhL4/EI.Therefore, substituting (4) into (1)–(3) will give(5)∂4w1∂x4+∂2w1∂t2+c∂w1∂t=α3Γ1∂2w1∂x2+α4VDC+VACcos⁡Ωt21-w12-α4VDC+VACcos⁡Ωt2d2/d1+w1-w22,(6)∂4w2∂x4+∂2w2∂t2+c∂w2∂t=α3Γ2∂2w2∂x2+α4VDC+VACcos⁡Ωt2d2/d1+w1-w22,(7)w10,t=0,w11,t=0,∂w1∂x0,t=0,∂w1∂x1,t=0,w20,t=0,w21,t=0,∂w2∂x0,t=0,∂w2∂x1,t=0,where the nondimensional parameters defined in (5)-(6) are given as follows:(8)α1=6Dh2;α2=6ε0L4Eh3D3;α3=6d1h2;α4=6ε0L4Eh3d13;c=c^L4TEI;Γ=∫01∂w∂x2dx;Γi=∫01∂wi∂x2dx.


## 3. Reduced-Order Modeling (ROM)

Reduced-order modeling (ROM) for the above described nonlinear problem can be obtained by discretizing (5) and (6) using the Galerkin expansion. Therefore, the deflection for the lower and upper microbeams will have the following forms, respectively:(9)w1x,t=∑i=1Npitϕix,(10)w2x,t=∑i=1Nqitϕix,where ϕi(x) are the linear undamped mode shapes of a clamped-clamped microbeam, which are orthogonal. After that, (5) will be multiplied by (1-w1)2(d2/d1-w1-w2)2 and (6) will be multiplied by (d2/d1-w1-w2)2 since this will help in reducing the computational cost. Next, the outcome will be multiplied by ϕj(x) and then integrated from x=0 to x=1 to give the following ROM equations:(11)∫x=01ϕjx1-∑i=1Npitϕix2∗d2d1+∑i=1Npitϕix-∑i=1Nqitϕix2∗∑i=1Npitϕiivx+∑i=1Np¨itϕix+c∑i=1Np˙itϕixdx=∫x=01ϕjx1-∑i=1Npitϕix2∗d2d1+∑i=1Npitϕix-∑i=1Nqitϕix2∗α3Γ1∑i=1Npitϕi′′x+d2d1+∑i=1Npitϕix-∑i=1Nqitϕix2α4VDC+VACcos⁡Ωt2-1-∑i=1Npitϕix2α4VDC+VACcos⁡Ωt2dx,∫x=01ϕjxd2d1+∑i=1Npitϕix-∑i=1Nqitϕix2∗∑i=1Nqitϕiivx+∑i=1Nq¨itϕix+c∑i=1Nq˙itϕixdx=∫x=01ϕjxd2d1+∑i=1Npitϕix-∑i=1Nqitϕix2∗α3Γ2∑i=1Nqitϕi′′x+α4VDC+VACcos⁡Ωt2dx.


## 4. Perturbation Analysis

Numerical solutions of the ROM dynamical equations provide good results for low nonlinearity values. However, when tuning the system’s geometrical parameters as well as the forcing amplitudes, one can increase the source of the nonlinearity, leading to nonlinear frequency responses which cannot be fully captured by ROM. As a result, perturbation theory was investigated to find the resonance of the double-microbeams especially for any high nonlinearity factors. The method of multiple scales was used by direct attack of the equations of motion [6, 24]. So, the variable for the time scales and its derivatives were defined as follows: (12)T0=t,D0=∂∂T0,T1=εt,D1=∂∂T1,T2=ε2t,D2=∂∂T2.Next, the damping coefficient c and the forcing amplitude VAC were scaled so that their nonlinearity effect will be balanced in the modulation equations [6, 24] as follows:(13)c=ε2c,VAC=ε3VAC,where ε is a bookkeeping parameter.We seek a solution for the lower and upper electrodes in the following form, respectively:(14)w1x,t,ε=ws1x+u1x,t=ws1x+εu11x,T0,T2+ε2u12x,T0,T2+ε3u13x,T0,T2+⋯,(15)w2x,t,ε=ws2x+u2x,t=ws2x+εu21x,T0,T2+ε2u22x,T0,T2+ε3u23x,T0,T2+⋯,where ws1 and ws2 are the static components of the microbeams’ deflection and u1 and u2 are their dynamic components for the lower and upper microbeams, respectively.Substituting (14) and (15) into (5) and (6) gives(16)∂4ws1∂x4+∂4u1∂x4+∂2u1∂t2+c∂u1∂t=α3Γws1+u1,ws1+u1∂2ws1∂x2+∂2u1∂x2+α4VDC+VACcos⁡Ωt21-ws1-u12-α4VDC+VACcos⁡Ωt2d2/d1+ws1+u1-ws2-u22,(17)∂4ws2∂x4+∂4u2∂x4+∂2u2∂t2+c∂u2∂t=α3Γws2+u2,ws2+u2∂2ws2∂x2+∂2u2∂x2+α4VDC+VACcos⁡Ωt2d2/d1+ws1+u1-ws2-u22.The electrical forces α4(VDC+VACcos⁡Ωt)2/(1-ws1-u1)2 and α4(VDC+VAC⁡cos(Ωt))2/(d2/d1+ws1+u1-ws2-u2)2 are expanded around (u1) and (u1-u2), respectively, using a Taylor-series expansion. Further, the expression of (VDC+VACcos⁡Ωt)2 will be approximated as VDC2+2VDCVACcos⁡Ωt. Using the aforementioned expressions and canceling the static equations and neglecting all small terms while equating the coefficients of like powers of ε, the following equations are obtained: For the upper microbeam, that is, (16),
(i)
orderε1:(18)L1u11,u21=0,(ii)
orderε2:(19)L1u12,u22=α3ws1′′Γu11,u11+2α3u11′′Γws1,u11+3α4VDC2u1121-ws14-3α4VDC2u112-2u11u21+u212d2/d1+ws1-ws24,(iii)
orderε3:(20)L1u13,u23=-2D0D2u11-cD0u11+2α3ws1′′Γu11,u12+2α3u11′′Γws1,u12+2α3u12′′Γws1,u11+α3Γu11,u11u11′′+6α4VDC2u11u121-ws14+4α4VDC2u1131-ws15+6α4VDC2-u11u12+u11u22+u12u21-u21u22d2/d1+ws1-ws24+4α4VDC2u113-3u112u21+3u11u212-u213d2/d1+ws1-ws25+2F~1cos⁡ΩT0-2F~2cos⁡ΩT0,whereL1 is a linear differential operator defined for any two functions (f and g) by(21)L1f,g=D02f+fiv-α3Γws1,ws1f′′-2α3ws1′′Γws1,f-2α4VDC2f1-ws13-2α4VDC2fd2/d1+ws1-ws23+2α4VDC2gd2/d1+ws1-ws23=0,and where(22)F~1=α4VDCVAC⁡1-ws12,F~2=α4VDCVAC⁡d2/d1+ws1-ws22.For the upper microbeam, that is, (17),
(i)
orderε1: (23)L2u11,u21=0,(ii)
orderε2:(24)L2u12,u22=α3ws2′′Γu21,u21+2α3u21′′Γws2,u21+3α4VDC2u112-2u11u21+u212d2/d1+ws1-ws24,(iii)
orderε3:(25)L2u13,u23=-2D0D2u21-cD0u21+2α3ws2′′Γu21,u22+2α3u21′′Γws2,u22+2α3u22′′Γws2,u21+α3Γu21,u21u21′′+6α4VDC2u11u12-u11u22-u12u21+u21u22d2/d1+ws1-ws24+4α4VDC2-u113+3u112u21-3u11u212+u213d2/d1+ws1-ws25+2F~2cosΩT0,whereL2 is a linear differential operator defined for any two functions (f and g) by(26)L2f,g=D02g+giv-α3Γws2,ws2g′′-2α3ws2′′Γws2,g+2α4VDC2fd2/d1+ws1-ws23-2α4VDC2gd2/d1+ws1-ws23.Assuming no internal resonances situations, the solutions of (18) and (23) are assumed to consist of only the directly excited modes, Φx, since the indirectly excited modes will die out in the presence of damping. As a result, the solutions of the dynamic components u11 and u21 of the lower and upper microbeams, respectively, are (27)u11x,T0,T2=AT2eiωT0+A-T2e-iωT0Φ1x,(28)u21x,T0,T2=AT2eiωT0+A-T2e-iωT0Φ2x,where A(T2) is a complex-valued function that is determined by imposing the solvability condition at the third order, the overbar denotes the complex conjugate, and ω and Φx are the natural frequency and corresponding eigenfunction of the directly excited modes, respectively. Here, the complex-valued function is considered to be the same for both microbeams, while the mode shapes will compensate for the resulting difference. This will help in obtaining the eigenvalue problem. Substituting (27) and (28) into (18) and (23) we obtain the following nonlinear coupled eigenvalue problem (EVP):(29)-ω2Φ1+Φ1iv-α3Γws1,ws1Φ1′′-2α3ws1′′Γws1,Φ1-2α4VDC2Φ11-ws13-2α4VDC2Φ1d2/d1+ws1-ws23+2α4VDC2Φ2d2/d1+ws1-ws23=0,-ω2Φ2+Φ2iv-α3Γws2,ws2Φ2′′-2α3ws2′′Γws2,Φ2-2α4VDC2Φ2d2/d1+ws1-ws23+2α4VDC2Φ1d2/d1+ws1-ws23=0.Next, the solutions u11 and u21 will be substituted in the second-order equations (ε2), which will give(30)L1u12,u22=A2e2iωT0+2AA-+A-2e-2iωT0h1x,L2u12,u22=A2e2iωT0+2AA-+A-2e-2iωT0h2x,where(31)h1x=α3ws1′′ΓΦ1,Φ1+2α3Φ1′′Γws1,Φ1+3α4VDC2Φ121-ws14-3α4VDC2Φ12d2/d1+ws1-ws24+6α4VDC2Φ1Φ2d2/d1+ws1-ws24-3α4VDC2Φ22d2/d1+ws1-ws24,h2x=α3ws2′′ΓΦ2,Φ2+2α3Φ2′′Γws2,Φ2+3α4VDC2Φ12d2/d1+ws1-ws24-6α4VDC2Φ1Φ2d2/d1+ws1-ws24+3α4VDC2Φ22d2/d1+ws1-ws24.Assume the particular solution of u12 and u22 as(32)u12=A2e2iωT0Ψ11+2AA-Ψ12+A-2e-2iωT0Ψ11,u22=A2e2iωT0Ψ21+2AA-Ψ22+A-2e-2iωT0Ψ21,where the functions Ψij are the solutions of the boundary value problems:(33)M1Ψ1i,2ωδ1i=h1x,M2Ψ2i,2ωδ1i=h2x,Ψ1j=0;Ψ1j′=0;Ψ2j=0;Ψ2j′=0;allevaluatedatx=0,x=1forj=1,2,and where δ is the Kronecker delta operator and the two linear differential operators M1 and M2 are defined as(34)M1Ψ1i,ω=Ψ1iiv-ω2Ψ1i-α3Γws1,ws1Ψ1i′′-2α3ws1′′Γws1,Ψ1i-2α4VDC2Ψ1i1-ws13-2α4VDC2Ψ1id2/d1+ws1-ws23+2α4VDC2Ψ2id2/d1+ws1-ws23,M2Ψ2i,ω=Ψ2iiv-ω2Ψ2i-α3Γws2,ws2Ψ2i′′-2α3ws2′′Γws2,Ψ2i+2α4VDC2Ψ1id2/d1+ws1-ws23-2α4VDC2Ψ2id2/d1+ws1-ws23.Introduce the detuning parameter σ to describe the nearness of the excitation frequency Ω to the fundamental natural frequency of the excited mode ω as(35)Ω=ω+ε2σ.The solutions of the first and second order are then substituted in (20) and (25), leading to(36)L1u13,u23=D02u13+u13iv-α3Γws1,ws1u13′′-2α3ws1′′Γws1,u13-2α4VDC2u131-ws13-2α4VDC2u13d2/d1+ws1-ws23+2α4VDC2u23d2/d1+ws1-ws23=-2iωA′Φ1x-icωΦ1xA+χ1xA2A-+F~1eiσ1T2-F~2eiσ2T2eiωT0+cc+NST;L2u13,u23=D02u23+u23iv-α3Γws2,ws2u23′′-2α3ws2′′Γws2,u23+2α4VDC2u13d2/d1+ws1-ws23-2α4VDC2u23d2/d1+ws1-ws23=-2iωA′Φ2x-icωAΦ2x+A2A-χ2x+F~2eiσ2T2eiωT0+cc+NST,where NST stands for the nonsecular terms and the functions χi are defined as (37)χ1x=2α3ws1′′ΓΦ1,Ψ11+2ΓΦ1,Ψ12+2α3Φ1′′Γws1,Ψ11+2Γws1,Ψ12+2α3Γws1,Φ1Ψ11′′+2Ψ12′′+3α3ΓΦ1,Φ1Φ1′′+6α4VDC21-ws14Ψ11Φ1+2Ψ12Φ1+12α4VDC21-ws15Φ13+6α4VDC2d2/d1+ws1-ws24-Φ1Ψ11-2Φ1Ψ12+Φ1Ψ21+2Φ1Ψ22+Φ2Ψ11+2Φ2Ψ12-Φ2Ψ21-2Φ2Ψ22+12α4VDC2d2/d1+ws1-ws25Φ13-3Φ12Φ2+3Φ1Φ22-Φ23;χ2x=2α3ws2′′ΓΦ2,Ψ21+2ΓΦ2,Ψ22+2α3Φ2′′Γws2,Ψ21+2Γws2,Ψ22+2α3Γws2,Φ2Ψ21′′+2Ψ22′′+3α3ΓΦ2,Φ2Φ2′′+6α4VDC2d2/d1+ws1-ws24Φ1Ψ11+2Φ1Ψ12-Φ1Ψ21-2Φ1Ψ22-Φ2Ψ11-2Φ2Ψ12+Φ2Ψ21+2Φ2Ψ22+12α4VDC2d2/d1+ws1-ws25-Φ13+3Φ12Φ2-3Φ1Φ22+Φ23.Hence, we need to eliminate the secular terms in (36) by seeking a particular solution free of secular terms in the form:(38)u13x,T0,T2=Υ1x,T2eiωT0,u23x,T0,T2=Υ2x,T2eiωT0.Then, (38) are substituted into (36). Next, the coefficients of eiωT0 in the resultant equations will be equated. The outcome of two equations will be then multiplied by two adjoints functions P(x) and Q(x), respectively, and then integrated by parts to transfer the derivatives from ϒ to the adjoints, leading to(39)∫x=01Υ1-ω2P+Piv-α3Γws1,ws1P′′-2α3ws1′′ΓP,ws1-2α4VDC2P1-ws13-2α4VDC2Pd2/d1+ws1-ws23dx+PΥ1′′′-P′Υ1′′+P′′Υ1′-P′′′Υ1-α3Γws1,ws1PΥ1′-P′Υ1+2α3ws1′P∫x=01ws1′′Υ1dx-2α3ws1′2PΥ1+2α3ws1′Υ1ΓP,ws1x=01+∫x=01Υ22α4VDC2Pd2/d1+ws1-ws23dx=∫x=01P-2iωA′Φ1x-icωΦ1xA+χ1xA2A-+F~1eiσ1T2-F~2eiσ2T2dx;∫x=01Υ2-ω2Q+Qiv-α3Γws2,ws2Q′′-2α3ws2′′ΓQ,ws2-2α4VDC2Qd2/d1+ws1-ws23dx+∫x=01Υ12α4VDC2Qd2/d1+ws1-ws23dx+QΥ2′′′-Q′Υ2′′+Q′′Υ2′-Q′′′Υ2-α3Γws2,ws2QΥ2′-Q′Υ2+2α3ws2′Q∫x=01ws2′′Υ2dx-2α3ws2′2QΥ2+2α3ws2′Υ2ΓQ,ws2x=01=∫x=01Q-2iωA′Φ2x-icωAΦ2x+A2A-χ2x+F~2eiσ2T2dx.By adding (39) and rearranging the terms, we can obtain the two adjoints equations governing both functions P(x) and Q(x) as follows:(40)-ω2P+Piv-α3Γws1,ws1P′′-2α3ws1′′ΓP,ws1-2α4VDC2P1-ws13-2α4VDC2Pd2/d1+ws1-ws23+2α4VDC2Qd2/d1+ws1-ws23=0;-ω2Q+Qiv-α3Γws2,ws2Q′′-2α3ws2′′ΓQ,ws2-2α4VDC2Qd2/d1+ws1-ws23+2α4VDC2Pd2/d1+ws1-ws23=0;which are the same as the first-order equations (29) so we call them self-adjoint equations. The solvability condition can be then obtained as(41)∫x=01-2iωA′PΦ1+QΦ2-icωAPΦ1+QΦ2+A2A-Pχ1+Qχ2+F~1Peiσ1T2+F~2eiσ2T2Q-Pdx=0.Using Euler’s formula for eiσ1T2 and eiσ2T2, then separating the real and imaginary parts of the solvability condition will give the following two coupled equations:(42)β′az1+z2+18ωa3S1+S2=-F1ωcos⁡γ1-1ωF22-F21cos⁡γ2,(43)a′z1+z2+c2az1+z2=F1ωsin⁡γ1+1ωF22-F21sin⁡γ2,where(44)z1=∫x=01PΦ1dx,z2=∫x=01QΦ2dx,S1=∫x=01Pχ1dx,S2=∫x=01Qχ2dx,F1=∫x=01F~1Pdx,F21=∫x=01F~2Pdx,F22=∫x=01F~2Qdx,γ=σT2-β.By squaring both sides of (42) and (43) and then adding the results, while remembering that (β′=σ-γ′), we get(45)a02σz1+z2+18ωa02S1+S22+c24z1+z22=F1+F22-F212ω2.Since the quality factor is related to the damping coefficient (c) by(46)Q=ωcz1+z2it can be replaced in the solvability equation (45) to examine the effect of the quality factor on the system dynamics.


## 5. ROM and Perturbation Results

To solve for the coupled equations of motion of both the lower and upper microbeams using the ROM, it is necessary to solve for the unknown functionspi(t) and qi(t) in (11) and then substituting them back into (9) and (10), which will give the deflection of both microbeams.For the method of multiple scales based perturbation technique, we start first by scaling the first-order coupled equations, (29) for both functions Ф1 and Ф2. Then, we substitute them into the boundary value problems equations (33) to solve for the unknown functions Ψ11, Ψ12, Ψ21, and Ψ22. Then, we solve for the unknown dynamic amplitude from (45) and the phase from either (42) or (43). Finally, the deflection of the lower and upper microbeams can be calculated using (14) and (15), respectively. More details about this numerical resolution technique can be found in [16]. Figure 2 is also accommodating the above numerical method description in summarizing the sequence of operations done in this particular method of multiple scales (MMS) technique.Figure 2Since each two consecutive natural frequencies of this system are very close to each other, the equations of the perturbation will be solved twice for each one of those two, and then the deflections around the two natural frequencies for each microbeam will be added to give the resultant deflections of both microbeams. A comparison of the obtained results for the two above described methods for the case of double-microbeams is shown in Figures3–5. The selected parameters are the same as assumed in the static analysis [18] and shown in Table 1, with the only difference being the applied voltages, VDC = 10 Volts and VAC = 0.5 Volts. In this present case, it can be seen from all figures that overall both assumed numerical methods are in good agreement for all different cases. While the perturbation analysis technique is less accurate than the ROM, it showed better global results in the presence of gross nonlinearity. It can predict all of the stable and unstable branches which the ROM is not capable of capturing. Also, it gives a better global and local dynamical picture about the microsystem’s behavior. Hence, the perturbation technique will be used in the following parametric study.Table 1
Assumed geometrical and material properties.
Parameter Value Parameter Value Beam length (L) 150µm Effective Young’s modulus (E) 166 GPa Beam thickness (h) 1.0µm Density (ρ) 2,332 kg/m3 Beam width (b) 4.0µm Air gap depths (d1 and d2, resp.) Cased1=d2 1.25µm 1.25µm Cased1<d2 1µm 1.5µm Cased1>d2 1.5µm 1µmFigure 3
Frequency response curves of the (a) upper microbeam and (b) lower microbeam assuming both ROM and perturbation method for the case ofd1=d2.
(a) (b)Figure 4
Frequency response curves of the (a) upper microbeam and (b) lower microbeam assuming both ROM and perturbation method for the case ofd1<d2.
(a) (b)Figure 5
Frequency response curves of the (a) upper microbeam and (b) lower microbeam assuming both ROM and perturbation method for the case ofd1>d2.
(a) (b)


## 6. Parametric Study

In order to investigate the effect of different physical and geometrical parameters on the dynamical amplitude of the double-microbeams based actuator, three different values were selected (low, medium, and high) for each chosen parameter; then the frequency response curve was generated for each microbeam to investigate the effects of each tuned parameter. The selected parameters are the length, the quality factor, and the DC and AC voltages. For all those cases, the depths of both air gaps are kept equal (d1=d2=1.25 μm) and the width (b) and height (h) for both microbeams are initially assumed to be 20 μm and 1.5 μm, respectively. The length (L) and the quality factor (Q) were equal to 210 μm and 50, respectively, for all the cases except when they are intentionally varied.


### 6.1. Effect of the Actuator Length (L)

Figure6 shows the amplitude versus the excitation frequency for three different lengths (210 μm, 410 μm, and 610 μm). From both figures, the significance of changing the length on the dynamical amplitude, even when small voltages were applied (i.e., VDC=2 Volts and VAC = 0.1 Volts), is clear. When the length is of a small value (L=210 μm), the two resonance peaks are very close to each other, which means that the applied voltages have a slight effect on the microstructure and so do not change its fundamental frequencies. However, when the value of the microbeam length is increased (L=410 and 610 μm), the distance between the two peaks becomes larger, meaning the fundamental frequency becomes lowered and the system is prone to experience a pull-in instability. Also, with a high value of the microbeam length, the nonlinearity becomes more dominant, producing a softening-type of behavior at the fundamental frequency for both microbeams.Figure 6
Effect of the microbeam length (L) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VDC=2 Volts and VAC=0.1 Volts.
(a) (b)


### 6.2. Effect of the Quality Factor (Q)

Figure7 shows the amplitude versus the excitation frequency for three different quality factors (50, 150, and 500) at high DC and AC voltages (i.e., VDC=10 Volts and VAC=0.5 Volts). It can be noticed that there is no effect for the quality factor on the resonance locations, which is reasonable since increasing the quality factor will decrease the damping which has no effect on the natural frequencies. However, the quality factor is shown to maximize the amplitude and with a high quality factor the nonlinearity becomes prominent (seehardening behavior in both figures).Figure 7
Effect of the quality factor (Q) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VDC=10 Volts and VAC=0.5 Volts.
(a) (b)


### 6.3. Effect of the Applied DC Static Load (VDC)

The effects of the applied DC voltage on the amplitude are shown in Figure8. Three different voltages were selected, which are 5, 10, and 15 Volts. At a low DC voltage, the fundamental frequency is not affected by the applied voltage, so it is very close to its consecutive one. However, an increase in the value of the DC voltage causes the fundamental frequency to shift from its initial position, and hence the system gets closer and closer to the pull-in instability. Since the term of the applied AC voltage is small there is no presence for the nonlinearity. So, we would like to investigate the nonlinearity effect for the applied DC voltage at a high AC voltage. As it can be shown from Figure 9, by increasing the value of the applied AC voltage to 5 Volts, the hardening-type nonlinearity appears at high values of the DC voltage. This indicates that the mid-plane stretching term is significant, since it is responsible for this nonlinearity.Figure 8
Effect of the applied DC voltage (VDC) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VAC = 0.1 Volts.
(a) (b)Figure 9
Effect of the applied DC voltage (VDC) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VAC = 5 Volts.
(a) (b)


### 6.4. Effect of the Applied AC Harmonic Load (VAC)

Figure10 shows the effect of the applied AC voltage on the amplitude. For this investigation, three different AC voltages were selected which are 5, 10, and 15 Volts. Unlike the DC voltage, the applied AC voltage has no effect on the natural frequencies, which is predictable since the dynamic loading has no effect on the natural frequencies. As shown in the figures, increasing the AC voltage increases the amplitude as well as producing a hardening type of nonlinearity.Figure 10
Effect of the applied AC voltage (VAC) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VDC = 2 Volts.
(a) (b)


## 6.1. Effect of the Actuator Length (L)

Figure6 shows the amplitude versus the excitation frequency for three different lengths (210 μm, 410 μm, and 610 μm). From both figures, the significance of changing the length on the dynamical amplitude, even when small voltages were applied (i.e., VDC=2 Volts and VAC = 0.1 Volts), is clear. When the length is of a small value (L=210 μm), the two resonance peaks are very close to each other, which means that the applied voltages have a slight effect on the microstructure and so do not change its fundamental frequencies. However, when the value of the microbeam length is increased (L=410 and 610 μm), the distance between the two peaks becomes larger, meaning the fundamental frequency becomes lowered and the system is prone to experience a pull-in instability. Also, with a high value of the microbeam length, the nonlinearity becomes more dominant, producing a softening-type of behavior at the fundamental frequency for both microbeams.Figure 6
Effect of the microbeam length (L) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VDC=2 Volts and VAC=0.1 Volts.
(a) (b)


## 6.2. Effect of the Quality Factor (Q)

Figure7 shows the amplitude versus the excitation frequency for three different quality factors (50, 150, and 500) at high DC and AC voltages (i.e., VDC=10 Volts and VAC=0.5 Volts). It can be noticed that there is no effect for the quality factor on the resonance locations, which is reasonable since increasing the quality factor will decrease the damping which has no effect on the natural frequencies. However, the quality factor is shown to maximize the amplitude and with a high quality factor the nonlinearity becomes prominent (seehardening behavior in both figures).Figure 7
Effect of the quality factor (Q) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VDC=10 Volts and VAC=0.5 Volts.
(a) (b)


## 6.3. Effect of the Applied DC Static Load (VDC)

The effects of the applied DC voltage on the amplitude are shown in Figure8. Three different voltages were selected, which are 5, 10, and 15 Volts. At a low DC voltage, the fundamental frequency is not affected by the applied voltage, so it is very close to its consecutive one. However, an increase in the value of the DC voltage causes the fundamental frequency to shift from its initial position, and hence the system gets closer and closer to the pull-in instability. Since the term of the applied AC voltage is small there is no presence for the nonlinearity. So, we would like to investigate the nonlinearity effect for the applied DC voltage at a high AC voltage. As it can be shown from Figure 9, by increasing the value of the applied AC voltage to 5 Volts, the hardening-type nonlinearity appears at high values of the DC voltage. This indicates that the mid-plane stretching term is significant, since it is responsible for this nonlinearity.Figure 8
Effect of the applied DC voltage (VDC) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VAC = 0.1 Volts.
(a) (b)Figure 9
Effect of the applied DC voltage (VDC) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VAC = 5 Volts.
(a) (b)


## 6.4. Effect of the Applied AC Harmonic Load (VAC)

Figure10 shows the effect of the applied AC voltage on the amplitude. For this investigation, three different AC voltages were selected which are 5, 10, and 15 Volts. Unlike the DC voltage, the applied AC voltage has no effect on the natural frequencies, which is predictable since the dynamic loading has no effect on the natural frequencies. As shown in the figures, increasing the AC voltage increases the amplitude as well as producing a hardening type of nonlinearity.Figure 10
Effect of the applied AC voltage (VAC) on the frequency responses of the (a) lower microbeam and (b) upper microbeam for VDC = 2 Volts.
(a) (b)


## 7. Conclusion

An investigation into the nonlinear dynamics of a multilayers based MEMS resonator made of clamped-clamped microbeams under electrostatic actuation was conducted. A nonlinear Euler-Bernoulli beam model was utilized accounting for both geometric (mid-plane stretching) and actuation nonlinearities. Two numerical methods were used to discretize the continuous beam model, which are perturbation analysis and Galerkin based reduced-order modeling (ROM). The results showed that the outcomes of the two methods are in good agreement. Although the ROM is more accurate in obtaining the behavior, it cannot capture all of the solution’s branches, where perturbation can be accomplished principally when the nonlinearity is becoming nontrivial.The influences of different geometrical and design parameters were examined. Based on the parametric study, the effects of the resonator parameters selection were shown to be significant on the MEMS multilayers resonator dynamical profiles. For example, increasing the resonator length (L) revealed a significant softening nonlinearity. The resonator quality factor (Q) decrease showed a huge effect in maximizing the resonator dynamical response. Regarding the electric force actuating loads, increasing the DC voltage (VDC) indicated a shift in the MEMS multilayers resonator resonant peaks away from each other, while producing a hardening-type nonlinearity for higher AC voltages. Furthermore, assume that higher AC load (VAC) produced more developed hardening behavior of the multilayers MEMS resonator. Consequently, interesting dynamical features can be achieved by controlling these resonator parameters, which may be useful in numerous MEMS resonating/sensing applications.



---

*Source: 1262650-2017-01-05.xml*