Abstract
We present a theoretical investigation of the dynamic behavior of a microelectromechanical system (in brief, MEMS) device modelled as a clamped-clamped microbeam subjected to electrostatic and electrodynamic actuation. We use the Galerkin projection technique to reduce the partial integro-differential equation governing the dynamics of the microbeam to a system of coupled ordinary differential equations which describe the interactions of the linear mode shapes of the microbeam. Analytical solutions are derived and their stability is studied for the simplest reduced-order model which takes into account only the first linear mode in the Galerkin procedure. We further investigate the influence of the first few higher modes on the Galerkin procedure, and hence its convergence, by analysing the boundaries between pull-in and pull-in-free vibrations domains in the space of actuation parameters. These are determined for the various multimode combinations using direct numerical time integration. Our results show that unsafe domains form V-like shapes for actuation frequencies close to the superharmonic, fundamental, and subharmonic resonances. They also reveal that the single first-mode reduced model usually considered underestimates the left branches and overestimates the right branches of these boundaries.
1. Introduction
Microelectromechanical systems (henceforth, MEMS) are structures that generally combine silicon based electrical and mechanical components at the scale of micrometers. MEMS technology emerged as a logical extension of microelectronics and integrated circuits technology [1–3] and was initially developed for realizing microsensors. Due to their small size, low fabrication cost, unique performance, and suitability for integration into complex functional engineering systems, they have become key components of many commercial systems. These include RF microresonators for wireless communication [4], accelerometers for airbag deployment in automobiles, ink jet printer heads, atomic force microscopes, optical switches, and chemical sensors [5]. The use of MEMS technology in the modelling of deformable mirrors for very high resolution giant telescopes has recently been addressed in adaptative optics [6]. The vibrating gyroscope is one of the MEMS devices that are commonly used for measuring angular velocity in area such as aviation, navigation, automotive, biomedicine, military affairs, and consumer electronics [7].
This large variety of MEMS can generally be classified according to their actuation mechanisms which can be pneumatic, thermal, piezoelectric, or electrical. The mathematical analysis of the latter actuation mechanism started in the late 1960s with the pioneering work of Nathanson and his coworkers [8–10]. They constructed and analyzed a model consisting of a parallel-plates capacitor of which one plate is fixed and the other is a movable mass attached to a spring. The fixed plate of the capacitor is connected with an alternating current (AC) voltage while the movable one is connected with a direct current (DC) polarization voltage. This model and its modified versions in which, essentially, the movable mass is replaced by a microbeam (either cantilever or fixed-fixed ends), have since then been studied in several aspects.
One of the most investigated problems in these studies is the phenomenon so-called pull-in instability. This phrase refers to the situation where the movable plate of the capacitor (microbeam or the mass) collapses on the fixed plate, leading to the failure of the MEMS devices. It is thus obvious that predicting the locations of its occurrence is of fundamental interest to MEMS designers. Studying the static (i.e., under the action of the polarization voltage only) deflection of the microbeam, Younis et al. [11] showed that the phenomenon corresponds to a saddle-node bifurcation. The dynamic pull-in phenomenon which accounts for the AC voltage in addition was reported and analyzed for switches actuated by a step voltage [12, 13] and various ramping rates [12]. Both studies [12, 13] indicate that the dynamic pull-in voltage can be as low as 91% of the static pull-in voltage. Nayfeh et al. [4] also investigated the pull-in instability phenomenon and identified three distinct mechanisms that lead to it. More recently, the influence of the size of the microstructure on the pull-in instability of a circular plate under hydrostatic and electrostatic actuation has been investigated by Mohammadi et al. [14]. They found that increasing the dimensionless internal length scale parameter or decreasing the applied hydrostatic pressures leads to higher values of the pull-in voltage.
The great majority of the investigations on driven MEMS and similar systems, as we could find in the literature, are performed under varied approximations. For instance, only a very low order truncated Taylor series expansion of the electric force is often effectively considered [15]. While considering the full nonlinear electric force, other works retain only the single first-mode model in their analysis [16–23]. The papers by Younis et al. [11] and Zhang and Zhao [24] are among the rare ones, to our knowledge, which considered both the full electric force and the effect of the higher modes. However their analysis treated only the static problem. We can thus see that the effect of higher modes on the full nonlinear dynamic behavior of MEMS under harmonic electrostatic forces has not been investigated. Yet, the dynamic of driven nonlinear systems is well known to be very sensitive to parameters and initial conditions, as illustrated by the large studies reported on nonlinear single-degree-of-freedom systems such as the Duffing oscillator, the van der Pol oscillator, or the simple pendulum. Investigations of driven MEMS dynamics that go beyond the suitability of the single first-mode reduced-order commonly assumed are therefore of interest. In this respect, we aim in this paper to study the dependence of the stability of MEMS devices modelled as clamped-clamped microbeams on their actuation parameters. Our approach is based on the direct numerical time integration following recent works [19, 22].
The paper is organized as follows. In the next section, we give the description of the model considered and the mathematical equation governing its dynamics. In Section 3, we use the Galerkin projection method to reduce the integro-differential equation of motion to systems of ordinary differential equations (ODEs). The simplest of these systems of ODEs is next solved analytically using the method of harmonic balance in Section 4. Then, in Section 5, we present the stability maps determined by numerical integration. In this numerical investigation we considered several modal combinations in order to appreciate the influence of higher order modes on the results. We end with our conclusions in Section 6.
2. Model and Equation of Motion
The device we consider in this work is a doubly-clamped microbeam illustrated schematically in Figure 1. It is modelled as a plate of length , width , and thickness submitted to a voltage , where , , and are, respectively, the DC polarization voltage and the amplitude and frequency of the applied AC voltage, and is the time. We assume that the microbeam is sufficiently narrow; that is, , so that its cylindrical deformations can be neglected. The transverse deflection of the microbeam is therefore assumed to be constant along its width and, thus, varies only along its length. Denoting the coordinate along the length of the microbeam by , the transverse deflection of the microbeam is governed by the equation [4, 11, 20] to which are associated the following boundary conditions: In (1a) is Young’s modulus, is the mass density, and are, respectively, the area and moment of inertia of the beam’s cross section, is the tensile or compressive axial load, is the capacitor gap width when , is the permitivity of free space, and is a viscous damping coefficient. Notice that the right-hand-side term of this equation represents an approximation of the parallel-plate electric force assuming complete overlap of the area of the microbeam and the stationary electrode. By introducing the dimensionless variables we reduce (1a) and (1b) to the following respective forms which are more convenient for our analysis:The time-rescaled voltage and the expressions of the new, dimensionless parameters , , , and in terms of the physical parameters of the system are given by
3. Galerkin Discretization
Following a standard approach [4, 5, 16, 20], we express the dimensionless transverse deflection of the microbeam as where with the th solution (in order of increasing magnitude) of the characteristic equation Kacem et al. [20] pointed out that for any integer , . It can be easily checked that the set of functions satisfy the orthonormalization conditions where is the Kronecker symbol. These functions as well as the characteristic equation (7) are obtained by seeking the solutions of (3a) and (3b) with in the form (5) and assuming , with . For the full nonlinear problem, the functions are generally not susceptible to such simple and exact analytical expressions. By multiplying (3a) by , then substituting (5) with (6) into the resulting equation, and finally integrating the outcome from to , we obtain a system of coupled second-order nonlinear ODEs governing the time evolution of the coefficients . Formula involving integrals and multiple summations for computing these ODEs can be found in [4, 20]. They are mostly suited for implementing the necessary operations in an algebraic manipulator such as MAPLE that we used for this work. Once these ODEs are derived, the coefficients of the expansions (3a) and (3b) can then be solved numerically or analytically using approximate methods from this system of nonlinear ODEs. In this paper we use consistently the method of harmonic balance to investigate analytically the single mode approximation in the next section and the numerical integration in Section 5 to analyze the multimode case.
4. Quasianalytical Approach for the First Modal Equation
Taking in (5) and applying the Galerkin procedure described above results to the following modal equation: The constant coefficients and introduced for convenience in this equations are given explicitly in Table 1. In spite of the least number of modes considered we obtain here a nonlinear ODE quite more complicated than most classically known ODEs, including the cubic Duffing oscillator, the cubic-quintic Duffing oscillator, the van der Pol oscillator, or the Helmhotz-Duffing oscillator. This implies some challenges in its investigation especially as the coefficients of the various nonlinear terms are all comparable in magnitude; see Table 1. So, for the case of a clamped-clamped beam symmetrically located between two electrodes, Rezazadeh et al. neglected all nonlinear terms in their single-mode reduction analysis which resulted in a standard linear Mathieu equation [15]. The variational iteration method [25] could then be used subsequently to determine the transition curves and solutions analytically.
In [20], Kacem et al. considered the alternative form of (9) where the whole equation is multiplied by the denominator of the last term. They invoked the method of averaging to analyze the slow modulation of the amplitude and phase of the fundamental response which they represented by an harmonic function. Such an approach will account only for the odd nonlinear terms of the equation. As it is well known [26, 27], a consistent analysis which can account for the quadratic terms , , , and requires the inclusion in the solution of at least a bias term and a second harmonic in addition to the fundamental one.
In this paper, we carry out an investigation beyond the linear regime while accounting for all nonlinear terms. To this end, we follow the approach used in [28]. Let us denote the stable equilibrium point of the unforced and undamped ODE associated with (9) by . This is the solution of the algebraic equation which further satisfies the condition We decompose the modal coefficient into a static components due to and a dynamic component due to . Thus, we write , where is a bookkeeping parameter introduced for the tractability of the calculations. We now substitute this into (9), which we next expand in power series of up to the third order. Finally we set and obtain a polynomial ODE of the form The various quantities are cumbersome functions of both and the coefficients , , and of (9). Their expressions are not presented here for the sake of conciseness. As pointed out by Rützel et al. [28] and Kacem et al. [20], (11) includes quadratic and cubic nonlinear terms, linear and nonlinear parametric excitation terms, and harmonic and superharmonic excitation terms. The resulting dynamic response can exhibit a vast range of nonlinear dynamic phenomena. However, the issues of existence and stability of periodic solutions of the microbeam are central to the operation of MEMS device and is focused on below.
To solve (11) we adopt one of the variants of the method of harmonic balance, namely, the so-called variable coefficients harmonic balance [29, 30]. According to this method, the sought periodic solution is expressed approximately as a truncated Fourier series with time-dependent coefficients. In permitting the amplitude of the periodic solution to vary with time, this method is able to determine the stability of the periodic solution. In a our case, we take We substitute the above in (11) and expand the outcome in Fourier series. From the coefficients of and , with we derive a system of five second-order autonomous ODEs for the coefficients of our approximate periodic solution. The stationary points of this latter system of ODEs correspond to periodic solutions of (11) and, thus, of the initial problem. The stability of these periodic solutions can also be assessed through the stability of the corresponding stationary points, which can be determined by analysing the eigenvalues of the jacobian matrix of the system. We refer the reader to [29, 30] for more details.
The whole procedure described above reduces the determination and characterization of the dependence of periodic solutions on the microbeam’s parameters to the study of fixed points of autonomous ODEs. This task is conveniently performed through computational tools such as Auto97 [31]. In this work we used a similar package called CL_MATCONT [32], which is developed to run within the widely used scientific computing environment MATLAB. Throughout this work we use the same numerical values for the microbeam’s parameters as in [33]: , , , and . The outcome of this study is summarized in Figures 2 and 3. Here the range of representation along the vertical axis corresponds approximately to the pull-in free domain. The solid lines denote stability and the dashed lines denote instability. Figure 2 depicts the response curves as a function of . It reveals a softening behavior characterized by curves which lean toward the left for the selected values . The responses are multivalued with two branches, each having both stable and unstable portions. The change of stability for each of the branches is a bifurcation whose nature and influence on the dynamics have been discussed with details in Nayfeh et al. [4]. The two branches move away from each other as is increased. This behavior is more acute for the main resonance. When is taken as control parameter, we obtain Figure 3. We can observe here a range of where there exists multistable response for frequencies relatively close to the resonance (up to three for ). However, this multistability disappears as the value of frequency gets more close to the resonant value.
(a)
(b)
5. Numerical Investigations
This section is devoted to the presentation of the stability maps delineating the pull-in from pull-in free vibrations in the plane. These maps are determined numerically using the following procedure. For each mode or combination of modes, the set of second-order differential equations obtained from the Galerkin projection is converted into a system of first order differential of the form where is a vector of components with the number of modes taken in the Galerkin projection and . Next the rectangular domain and is subdivided into mesh points with a resolution of . For each mesh point, we use the Runge-kutta of order 4 as implemented in the ode45 routine of the MATLAB package to solve (13). The initial conditions are chosen as the physically stable equilibrium position (static deflection) of the microbeam and are solved for numerically from (13) with , using the Newton-Raphson algorithm. Finally the mesh point is deemed to correspond to pull-in free vibration if the time-integration of the ODEs system can be performed for a time length of 500 periods of the driving voltage without observing the pull-in condition described below. This type of analysis can be performed also in the state space of the system. This is the case of [19, 22] where basin of attraction and their erosion are used to characterize the integrity of the system studied.
For the single, first-mode approximation, the modal equation (9) can be seen as that of a single degree-of-freedom particle embedded in the potential depicted in Figure 4 for our choice of the parameters’ values. This potential has two mathematically stable points at and and one unstable point at . Around the stable point , the microbeam deflection can grow up to one. In fact the microbeam deflection for the first mode is . A simple graph of shows that is maximum value is . This implies that, for the microbeam deflection to be less than one, must be less than for all . Motions of the particle either with initial conditions chosen in the right wells of the potential or which lead to its escape from left wells cannot obviously fulfil this constraint. Therefore the stability condition is that the particle is confined in the potential wells that contain the stable position which is close to zero. This definition can be appropriately extended to the case where several modes are included in the Galerkin projection.
Figure 5 presents a typical stability map for our MEMS. The blackened regions are pull-in free while the white generally -shaped areas are those leading to pull-in instability. The largest region of instability is near the system’s natural frequency as given in [33] for the same MEMS. This unstable region characterizes the pull-in phenomenon at the primary resonance. In the same figure, we have other unstable regions at , and . The first two correspond to superharmonic resonances of respective orders and while the last corresponds to a subharmonic resonance of order 2.
We provide examples of the time history and phase diagram of the modal coefficient in both the pull-in free and pull-in situation in Figures 6 and 7. They are plotted, respectively, for points labelled and in Figure 5. We see that at , the time evolution maintains so that the microbeam deflection fulfils . For this particular pair of parameters’ values the time evolution is multiperiodic with some sort of modulation. For point , Figure 7(a) shows on the contrary that evolves chaotically with time. We emphasize however that it is the fact that for some times, which implies that one would have , which determines this point as unstable.
(a)
(b)
(a)
(b)
We now turn to the investigation of the influence of higher modes on the Galerkin projection. To this end, we repeated the procedure presented above for several combinations of the modes 1, 2, 3, 4, and 5. The results generally look like the one of Figure 5. To ease the comparison of these results, we extracted the boundaries between the pull-in and pull-in free regions for the various combinations considered and superimposed them on the same graphs.
We consider first the second linear mode by taking in (9). Then the Galerkin technique yields the following system of ODEs: whereWe display the corresponding boundaries in Figure 8. It appears in this case that the two boundaries are strictly identical. This indicates that second mode does not contribute to the dynamics. In their study, Younis et al. considered only the odd modes [11]. Zhang and Zhao [24] reported that even modes do not largely influence the static pull-in instability. We should keep in mind, however, that we have used here a specific set of initial conditions associated with the static deflection of the microbeam. The situation might not be the same for different choices of the initial conditions.
When the third-order linear mode is included in addition of the first two, we obtain Figure 9. As one can see, the effect of the third mode is to shift the basin of stability to the right and slightly increase the instability regions by lowering its threshold around the resonances. Here also, by performing the analysis with only the first two odd modes and , we obtained exactly the same figure. This indicates that the even mode does not contribute to the dynamic of the MEMS device when it is operated from its static deflection. In fact, we noticed in our numerical investigations that this is also the case for any combination of modes 1–5 which includes the second even mode . Henceforth, our labelling of the various curves displayed, especially in Figures 10(a)–10(c), accounts only of the modes which have been found to influence the boundaries for our choice of initial conditions.
(a)
(b)
(c)
At last, we show in Figures 10(a)–10(c) the modifications induced by the inclusion of the third even mode . It is observed that one influence of the higher modes is in general to shift the instability regions toward the right. It also appears from these figures that the single first-mode approximation gives the lowest estimation of the left branches and the highest estimation of right branches of the -shaped bifurcation curves. These figures seem to indicate a monotonous ordering of the transition curves with the number of modes taken in the Galerkin projection only for the left branches. They also suggest that the right branches are bracketed by the estimates of the first-mode and the first two odd-mode approximations.
6. Conclusion
In this paper we have been concerned by the dynamics of a MEMS device governed by an integro-differential equation. We used the Galerkin projection technique to transform this integro-differential equation into various systems of ODEs. In the case of the single first-mode reduction, the response curves and their stability have been obtained as a function of the alternating voltage based on the harmonic balance method. Through the numerical integration of the modal ODEs, we have performed a systematic investigation of the dependence of safe dynamics of the MEMS device on the parameters of the actuation, using various combinations of modal approximation. Our study reveals that unsafe domains form -like shapes for actuation frequencies close to the superharmonic, fundamental, and subharmonic resonances. The fundamental natural frequency of a resonant microbeam is the most sensitive to electrodynamic excitation parameters (frequency and amplitude). The various modal combinations considered in our work do not show a clear monotonous ordering of the transition curves with the number of modes taken in the Galerkin projection.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors acknowledge the critical reading of the paper by Professor Anindya Chatterjee of the Department of Mechanical Engineering, India Institute of Technology, Kanpur 208016, India.