Nonlinear Behavior of Electrostatically Actuated Microbeams with Coupled Longitudinal–Transversal Vibration

Microelectromechanical switch has become an essential component in a wide variety of applications, ranging from biomechanics and aerospace engineering to consumer electronics. Electrostatically actuated microbeams and microplates are chief parts of many MEMS instruments. In this study, the nonlinear characteristics of coupled longitudinal–transversal vibration are analyzed, while an electrostatically actuated microbeam is designed considering that the frequency ratio is two to one between the first longitudinal vibration and transversal vibration. The nonlinear governing equations are truncated into a set of coupled ordinary differential equations by the Galerkin method. Then the equations are solved using the multiple-scales method and the nonlinear dynamics of the internal resonance is investigated. The influence of bias voltage, longitudinal excitation and frequency detuning parameters are mainly analyzed. Results show that using the pseudo-arclength continuation method, the nonlinear amplitude–response curves can be plotted continuously. The saturation and jump phenomena are greatly affected by the bias voltage and the detuning frequency. Beyond the critical excitation amplitude, the response energy will transfer from the longitudinal motion to the transversal motion, even the excitation is employed on the longitudinal direction. The large-amplitude jump of the low-order vibration mode can be used to detect the variation of the conditions or parameters, which shows great potential in improving precision of MEMS switches.

Electrostatically actuated microbeams and plates exhibit significant nonlinearities due to the action of the electrostatic force, as the electrostatic force is inversely proportional to the square of the distance between the electrode and the host structure. On account of nonlinearities, a collapse of the movable structure occurs at a critical voltage (pull-in instability), and the phenomenon can be used as change of ON or OFF state [18][19][20][21][22][23][24][25]. Younis [26,27] used analytical approaches to investigate the behavior of electrically actuated microbeam-based MEMS, and proposed a reduced-order model in nonlinear dynamic analysis of MEMS. Alsaleem [28] showed experimental and theoretical investigations of dynamic pull-in of electrostatically actuated resonators, and the influences of different conditions were investigated. Ouakad and Younis [29,30] studied the dynamic behavior of clamped-clamped micromachined arches with a Galerkin-based reduced-order model, a variety of nonlinear phenomena was analyzed in detail, such as hysteresis, softening behavior, dynamic snap-through, and dynamic pull-in. Wang [31] presented a size-dependent model for electrostatically actuated microbeam-based MEMS using strain gradient elasticity theory, in the study the relation of the geometry size and the normalized pull-in voltage were analyzed in detail.
Early studies mainly focused on the static and dynamic behavior of microbeam-based MEMS. Li [32] gave monostable dynamic analysis of microbeam-based resonators using an improved reduced-order model, results showed that the ratio of the gap width to microbeam thickness affected large-amplitude vibration significantly. Ghayesh [33] investigated the nonlinear size-dependent behavior of an electrically actuated MEMS resonator based on the modified couple stress theory and a high-dimensional reduced-order model. Based on reduced-order model, Caruntu [34] researched parametric resonance of microelectromechanical (MEMS) cantilever resonators under soft damping and soft alternating current (AC) electrostatic actuation. Kacem [35] constructed a comprehensive multiphysics model for electrostatically actuated clamped-clamped resonators based on the Galerkin decomposition method coupled with the averaging method. Feng [36] showed detailed analysis of class of bipolar electrostatically actuated micro-resonators, and the electrostatic force nonlinearity, neutral surface tension, and neutral surface bending were considered. Han [37] investigated the static and dynamic characteristics of a doubly clamped micro-beam-based resonator driven by two electrodes, nonlinear dynamic analysis for two cases were analyzed, including that the origin of the system was at a stable center or an unstable saddle point.
To get more precise results, Muldavin [38] presented an accurate model of the switching mechanism of MEMS switches based on an electro-mechanical analysis, varying force and damping versus position (time) were taken into consideration. Ghayesh [39] investigated the size-dependent dynamical performance of a microgyroscope using the modified couple stress theory, and size-dependent frequency-response curves for the sense and drive directions were obtained to analyze the dynamical performance of the microgyroscope. Li [40] used the modified couple stress-based strain gradient theory to construct a unified nonlinear model for an electrostatic MEMS microbeam capacitive switch, the quasistatic and dynamic behavior of which were studied systematically. Leus [41] applied energy methods to study the undamped dynamic response of electrostatic MEMS switches under a step-function voltage. Larose [42] analyzed the influences of inertial effects, structural, air damping (squeeze-film damping) and the impact behavior of the microbeam. Lin [43] mainly investigated the action of Casimir effect on the pull-in parameters of nanometer switches. Shear deformation was taken into account in the analysis of microbeams with geometric nonlinearities [44][45][46][47][48][49].
The dynamical behavior of microelectromechanical systems (MEMS) plays a pivotal role in determining their performance, such as natural frequency [50][51][52][53], critical pull-in voltage [54][55][56], bifurcations [57,58]. Dynamical bifurcations have been taken into account in MEMS sensors for mass detection [59][60][61], which may lead to new mass measurement architectures to miniaturized mass spectrometers [62]. Studies on coupled vibrations have introduced a host of nonlinear phenomena in micro-and nano-scale resonators for their use to reveal the mechanism of the complex dynamic behaviors. Alkharabsheh [63] investigated the effect of axial forces on the static behavior and the fundamental natural frequency of electrostatically actuated MEMS arches. Li [64] investigated the nonlinear characteristics of microbeam-based resonators under higher-order modes excitation.
It can be concluded that coupled vibration behaviors caused by internal resonance are gradually considered, and the jump phenomenon is also focused on [65][66][67]. The phenomenon is interesting, but the application of this phenomenon is scarce. Classical MEMS switches are usually designed based on pull-in instability, here this study aims to investigate application of jump phenomena in nonlinear bifurcation in improving the precision of MEMS switches. A nonlinear model for a microbeam is established, taking electrostatic force, axial load, coupled longitudinal-transversal vibrations and mid-plane stretching into consideration. A single mode approximation is derived and the physical condition for the jump is determined, then the influence of longitudinal excitation and frequency detuning parameters are analyzed qualitatively. Pseudo-arclength continuation method and a direct time-integration technique are used to investigate the nonlinear dynamic behavior of the electrostatically actuated microbeam due to the axial load. Using numerical method, frequency-amplitude responses and force-amplitude responses are given to demonstrate the nonlinear characteristics.
The rest of the paper is arranged as follows: In Section 2, mathematical model is given using Galerkin's method, and the multiple-scales method is used to derive an approximate average equation. Static analysis with convergence analysis and perturbation analysis are shown in Sections 3 and 4, respectively. A brief introduction of pseudo-arclength continuation is given in Section 5. Section 6 gives resonant conditions for this kind of electrostatically actuated microbeam. In Section 7, nonlinear dynamical properties of the resonant microbeam are investigated in detail. Several conclusions are summarized in the last section.

Mathematical Model
A slender microbeam under the action of electric field and axial load is considered, as shown in Figure 1. The total length of the beam is L, with a thickness of h and width of b. The distance between the fixed electrode and the microbeam is d. On the right end of the beam, an axial force is applied. Two electric field V dc is applied to detune the stiffness of the switch. In the following derivation, same assumptions from Reference [33] are to be applied:(1) the Euler-Bernoulli beam theory is applied without regard to shear deformation effect and rotary inertia; (2) the microbeam and the parallel-plate electrodes have a complete overlapping area; (3) the geometric nonlinearity is caused by the mid-plane stretching of the microbeam; (4) the transversal deflection is constant along the width. By Hamilton's principle and nonlinear Euler-Bernoulli beam theory, the coupled longitudinal-transversal vibration of the microbeam can be written as [32,52]: whereū is the longitudinal displacement andw is the transversal displacement; E and I z are the Young's modulus and moment of inertia of the cross section of the microbeam, ρ is the material density and A is the area of the cross section. An axial forceF cos(Ω 1 t) is applied at the right end of the microbeam. The terms on the right-hand side of Equation (2) represent the axial force, mid-plane stretching effects, and two parallel-plate electric actuation, respectively. In the approximated electric forces, C is the fringing field coefficient and 0 is the dielectric constant of the gap medium. It should be mentioned that the fringing field always exists in the electrostatically actuated microbeams. However, in this study, we aim to investigate the bifurcation phenomena and its potential in improving the precision of MEMS switches. Considering that the response trend will not be influenced by the fringing effect, the fringing effect is neglected for the sake of simplicity, and here C = 1. The boundary conditions for this microbeam are: To simplify the dynamic equations, the following dimensionless variables are applied: EI z is the time scale. Substituting nondimensional variables into Equations (1) and (2) yields where the new parameters are given as: To analyze the dynamic equations qualitatively, a reduced-order model is firstly derived using Galerkin's method, transforming the complex equations into a multi-degree of freedom system. The solutions of the equations can be presented by mode shapes and generalized coordinates: where φ k (x) and ψ k (x) are mode shapes for transversal vibration and longitudinal vibration, respectively; q k (t) and p k (t) are generalized coordinates for transversal vibration and longitudinal vibration, respectively. The linear undamped mode shape of the transversal vibration of the straight beam is: where λ k = β 2 k EI z ρA , and the following relationship should also be satisfied for β k : For the transversal vibration of the straight beam, one characteristic of the linear undamped mode shape is: For longitudinal vibration of a fixed-free beam, the linear undamped mode shape has the form: where α k = (2k+1)π 2L E ρ and β k L = (4.7300, 8.8532, 10.9956, · · ·). In Galerkin's method, the linear mode shapes of the beam can be used as basis functions, and the differential Equation (5) have the following equivalent form: For small w, the electric actuation force is simplified using Taylor approximation, and carried out as follows [68]: Using Equations (11) and (14), multiplying Equation (6) by φ, and integrating the outcome from x = 0 to 1, yields: Taking the orthogonality of mode shapes and Equation (11) into consideration, and only use the first order mode function, the coupling transversal-longitudinal vibration equations are obtained: where

Static Analysis and Convergence Analysis
An example of eigenfrequency analysis in the steady state are given at first. Here, results by single mode Galerkin approximation, finite element simulations, and results from previous literature are performed. The parameters from previous literature is applied in the calculation, and the calculated simulation results are shown in Figure 2, one can see that the results agree well with each other. To activate the internal resonance between two vibrational modes, dimensions of the switch need carefully designed. A two-to-one frequency ratio between the longitudinal mode and transversal mode can cause the energy exchange between the two modes. However, a perfectly 2:1 frequency ratio is not always effective and feasible in practical engineering, as small changes in structural dimensions and material properties can deviate the nominal 2:1 ratio. In the numerical simulation, the material and geometric parameters for the microbeam are designed as E = 169 GPa, ρ = 2300 kg/m 3 , L = 24 µm, h = 2.92 µm, g = 0.5 µm and b = 10 µm. Using finite element software ANSYS R15.0, modal analysis is carried out, and the first transversal mode shape and longitudinal mode shape are plotted in Figure 3. The natural frequency of the first longitudinal mode is 8.93 × 10 7 Hz and the natural frequency of the first transversal mode is 4.47 × 10 7 Hz, therefore a frequency ratio near two to one is obtained. Then a bias DC voltage between the electrodes and microbeam is applied, resulting in two electrostatic forces which will change the mechanical stiffness of the microbeam and consequently a frequency ratio of 2:1 is ensured. To verify the validity of the coupled Equations (16) and (17), the convergence analysis is carried out, comparing the results obtained by the reduced-order model and the results by direct numerical simulation of Equations (5) and (6). The direct numerical simulation is given by differential quadrature method [32]. The force-amplitude curves are simulated as shown in Figure 4, and three cases with different bias voltages are calculated. The results predicted by those two methods agree very well for small deformations. Therefore, the single mode model can describe the dynamic behavior of the system with small deformation assumption. It is interesting that two kinds of force-amplitude responses are obtained, one for V dc = 30 V with a jump in the amplitude response, and another for V dc = 60 V and V dc = 90 V with continually increasing amplitude, and the difference will be explained in Section 7. Amplitude of transversal vibraiton Single mode DQM Figure 4. Comparison of the force-amplitude response curves obtained by single mode Galerkin approximation and differential quadrature method.

Perturbation Analysis
In this section, the method of multiple-scales is to be applied to analyze the steady-state responses of the microbeam around equilibrium position. Before the analysis, a small nondimensional bookkeeping parameter is introduced to indicate the significance of each terms in Equations (16) and (17). Taking the electrostatic force term F 1 = O( 2 ) into consideration, the coupling vibration equations become:p q + 3 µ 1q + λ 2 q + γ 2 pq + γ 3 q + γ 4 q 3 = 0 From Equation (19), one can find that there is no transversal vibration when the amplitude of external force F 1 is small or driving frequency Ω is far away from two times of the first nature frequency of the transversal vibration. For the two-to-one internal resonance, two parameters ∆ and δ are introduced to describe the nearby vicinity of the resonance.
Based on multiple-scales method, the approximate solutions of the longitudinal vibration and transversal vibration can be defined in the following forms: where T n = n t represents different timescale. Substituting the approximate solutions into the dynamic equations and equating the coefficients of equal powers of , a series of second-order ordinary differential equations will be obtained.
The general solutions of Equation (24) can be written as whereĀ andB are complex conjugate functions of A and B, and j represents imaginary component. For convenience, the parameters in above solutions can be represented in polar form: where a 1 and b 1 represent the amplitudes of the first longitudinal vibration and the first transversal vibration, respectively. Substituting Equations (29) and (30) into Equations (25)-(28), the following secular terms are yielded: The papameters ϕ, χ, k 1 , k 2 and k 3 are defined as: To obtain the steady-state solutions, one can set all the time derivatives (ȧ 1 ,φ,ḃ 1 ,χ) to be zero and all the time-dependent variables to be constants in Equations (31)-(34), then solve the algebraic equations with the right-hand side zero. The stability of steady-state solutions can be determined by the eigenvalues of Jacobian matrix of Equations (31)-(34) near the steady-state solutions. If all the eigenvalues of the Jacobian matrix have non-positive real parts, the system is stable otherwise the system is unstable. Finally, the steady-state frequency responses are determined by the following frequency-response equations: where s = δ + ∆ + k 2 a 2 1 + k 1 b 2 1 . With the Routh-Hurwitz criterion, stability of the system can be easily analyzed by the characteristic polynomial. Routh-Hurwitz stability criterion supplies a set of necessary and sufficient conditions for the accurate delineation of the relevant parameter space into stable and unstable regions. To determine the stability of the periodic solution, the Jacobian matrix of Equations (31)-(34) is obtained as shown in Equations (31)- (41), then the stability of the solutions is analyzed by the Routh-Hurwitz criterion.

Pseudo-Arclength Continuation
To analyze the steady states, Equations (39) and (40) are calculated by a two-step pseudo-arclength continuation method. Firstly, one advances along a branch of steady states with a varied parameter, then the linear stability analysis of the most recently computed steady state is carried out [69,70].
Equations (39) and (40) can be written in the matrix form: where u = (a 1 , b 1 ) and ς is the variation parameter. Function Φ(u, p) contains two functions here: The pseudo-arclength continuation is naturally a predictor-corrector method, parametrizing branches of solutions Γ(s) = (u(s), ς(s)) with an arclength parameter s. For every given solution (u 0 , ς 0 ) and the next solution (u, ς), the following relationship should be satisfied for a small step size ∆s: Furthermore, a normalization equation is supplemented to solve the extended system [70]: |Γ 0 (s)|= 1 (46) whereΓ 0 = (u 0 ,ς 0 ) represents the normalized direction vector of the solution family Γ(s) at (u 0 , ς 0 ). To compute the normalized direction vector Γ 0 , one can solve the following equation: where Φ u is the derivative of Jacobian matrix to parameter u and Φ ς is the derivative of Jacobian matrix to parameter ς at point (u 0 , ς 0 ). The predictor solution of is given by In corrector algorithm, the predictor solution (u 0 , ς 0 ) is projected back to the branch in a direction orthogonal to the tangentΓ 0 . The Newton-Raphson iterations is used here, and the iteration is performed as: where r k = ∆s −u T 0 (u k − u 0 ) −ς 0 (ς k − ς 0 ). After getting the new iteration direction, (u k , ς k ) is updated by (u k+1 , ς k+1 ).
In practical calculation, it is sometimes better to solve two n × n linear systems instead of directly solving, namely Finally, the new iteration direction can be obtained using the new parameters z 1 and z 2 : The advantage of pseudo-arclength method is that the Jacobian matrix of the extended system has rank n [69], even at folds where Φ u becomes singular. Hence, using pseudo-arclength method, one can plot the frequency responses curves easily without piecewise process, and the stable and unstable solutions around folds can be computed continuously.

The Resonant Condition
The problem of Equations (1) and (2) can be classified as a parametrically excited system. The resonance occurs if the excitation frequency of the longitudinal vibration approaches twice of the any frequency of the transversal vibration, and the resonance curve typically displays a trivial solution [58]. To determine the critical state of the dynamic system, it is advantageous to transform the general solutions from polar coordinates b 1 and ϕ to rectangular coordinates η and : Substituting Equation (56) into Equations (31) and (32), the derivatives of η and can be obtained: Similar to above introductions, the stability of steady-state solutions of Equations (57) and (58) can also be determined by the eigenvalues of the linearized coefficients matrix (i.e., Jacobian matrix) near the steady-state solutions. If all the eigenvalues of the Jacobian matrix all have negative real part at an equilibrium point, the point is asymptotically stable, the system is stable, otherwise at least one eigenvalue has positive real part, the equilibrium is unstable. The linearized coefficients matrix has the following form: According to Equation (59), the local stability can be analyzed using the trace and determinant of the Jacobian matrix. As c n is always larger than zero here, the critical point can be found while the determinant of the Jacobian matrix is zero. With Det(J) = 0, the threshold for a 1c is obtained: The threshold implies that: if a 1 > a 1c , the transversal vibration may occur, else if a 1 < a 1c there is no transversal vibration. For an a 1c in physical meaning, the value in the radical sign must be non-negative. Therefore, the physical condition for the transversal vibration is: The basic physical condition is also the critical condition for modal coupling vibration. As the amplitude of the longitudinal vibration increases exceeding the critical value a 1c , the energy transfers from the longitudinal mode to the transversal mode. The stability analysis of the nontrivial solution should be given as the nontrivial solution branches bifurcate. For this parametrically excited system, Near the critical points, the supercritical Hopf bifurcation leads to stable branches, while subcritical Hopf bifurcation leads to unstable branches. To characterize the stability of periodic vibration, Hopf bifurcation of critical points is to be studied. Substituting the critical solution a 1c into Equation (39), the following discriminant can be obtained: If Λ < 0, the subcritical Hopf bifurcation occurs, and the jump phenomenon will appear in the transversal mode when the amplitude of the longitudinal vibration exceeds the critical value. On the contrary, if Λ > 0, the supercritical Hopf bifurcation occurs, and the longitudinal vibration only induces the small transversal vibration, even the amplitude of the longitudinal vibration is large. For Λ = 0, the threshold of amplitude of the longitudinal vibration is minimum, meaning that a relatively smaller electrostatic force could motivate the transversal vibration.

Dynamic Analysis
To further research the nonlinear dynamical behavior under different Hopf bifurcation parameter range, the influences of electrostatic force and detuning frequency on the system are introduced. In this section, we study the complex dynamical behaviors of electrostatically actuated microbeam using pseudo-arclength continuation method and some interesting phenomena are obtained. The dynamic responses of the system are simulated by Crank-Nicolson method, which is numerically stable.
Firstly, the bifurcation behavior of the longitudinal-transversal coupling vibration is investigated. Figure 5 displays variation of the bifurcation behavior versus V dc and δ. It is found that a larger DC voltage corresponds to a larger δ, indicating that with a larger DC voltage, the modal coupling coefficient is reinforced, and nonlinear modal interactions is enhanced. Here we can explain the phenomenon in Figure 4: in the simulation δ = 0.1 is used, with V dc = 30 V, subcritical bifurcation occurs while V dc = 60 V and V dc = 90 V supercritical bifurcation occurs, hence two kinds of force-amplitude responses are obtained. In Figure 6, V dc = 140 V and δ = −0.1 is employed, and subcritical Hopf bifurcation occurs. Hence as F 1 increases, the jump phenomenon appears in the second-order mode and the amplitude of the second mode is much larger than that of the third mode. The phenomenon of saturation can be found for the longitudinal vibration, and the energy transfers from longitudinal mode to transversal mode. In Figure 6, as the external force increases from zero, the amplitude of the longitudinal vibration increases along the natural relationship A = , and the amplitude of the transversal vibration still keeps zeros. It is interesting that beyond a critical value of F 1 , the solution of longitudinal vibration loss stability, and another branch of solution dominates the motion. Then the amplitude of the longitudinal vibration becomes independent of the amplitude of the external force, i.e., the saturation occurs. It is interesting that the jump phenomenon occurs in transversal mode, and the jump amplitude is very large. Therefore, taking advantage of the nonlinear features of the frequency response, a large-amplitude jump will provide an effective approach to improve the precision of MEMS switches. In Figure 7, supercritical Hopf bifurcation occurs while V dc = 140 V and δ = 0.1 is adopted. From Figure 7, the phenomenon of saturation of longitudinal mode can also be found. However different from the plots in Figure 6, as the external force increases from zero, the amplitude of the longitudinal vibration increases as well, while the amplitude of the transversal vibration keeps zeros. Beyond a critical value of F 1 , the amplitude of transversal vibration experiences a drastic rise first, and then increases in a slow rate continually.  Figure 8 shows the amplitude responses of the coupled vibration as function of internal resonance detuning parameter δ. A constant bias voltage V dc = 140 V is applied to supply an electrostatic force. The frequency sweep response shows a slight asymmetry. Before reaching the threshold, the longitudinal vibration operates in its linear regime and the amplitude of longitudinal vibration increases linearly, while there is no transversal vibration. Once reaching the certain threshold, the solution of longitudinal vibration turns to another branch of solution. Meanwhile an amplitude jump of transversal vibration occurs, and the vibration energy transfer from the longitudinal mode into the transversal mode. In the design of switches, the jump phenomenon can be used to change the ON and OFF state, as the amplitude jump is very large, the precision and switch speed can be greatly promoted.
As introduced above, the jump phenomenon is used and the parametric excitation vibration only in the subcritical Hopf regime is taken into consideration in this section. The parametric force-response curves obtained by pseudo-arclength method for different values of DC voltages are plotted in Figure 9. In the calculation, δ = −0.15 is used and V dc ranges from 50 V to 200 V with an interval of 50 V. It is found that with a larger DC voltage, a larger excitation F 1 is needed to realize the jump phenomenon of the transversal mode. As also seen from the plots, with a larger DC voltage, the amplitude of the jump is larger as well, demonstrating that the performance of the switch can be improved by increasing the DC voltage. The force-response curves for different detuning parameter δ are plotted in Figure 10. Here V dc = 150 and δ = (−0.05, −0.1, −0.15) is used. The plots are similar with that in Figure 10. However, the point of the limit point bifurcation and the period-doubling bifurcation are all delayed while δ is smaller, and with a smaller δ the amplitude of the jump is larger. Both in the two figures, a larger saturation value can be found with a larger jump, illustrating that the energy is transferred between the longitudinal-transversal modes interaction.  Figure 11 shows the frequency sweep response of transversal vibration for different V dc . The values of V dc is defined as 50 V, 100 V, 150 V, and 200 V. The whole frequency sweep response curves are similar with the plots in Figure 8. It is visualized that the unsymmetrical configuration can still be found, and the asymmetry is strengthened by a larger DC voltage. Here only the curves for transversal mode are given. The plots show that with a lager bias voltage, the left jump is also bigger while the right jump is smaller. Both the left and right threshold shift right while the bias voltage increases.

Conclusions
To extend the application of bifurcations of the coupling longitudinal-transversal vibration in the design of MEMS switches, this paper mainly focused on the nonlinear characteristics of an electrostatically actuated microbeam, and the coupled longitudinal-transversal vibration was analyzed in detail. A reduced-order model was obtained by Euler-Bernoulli beam theory and Galerkin's method, and the obtained single mode Galerkin model was verified by differential quadrature method, demonstrating that the single mode model could describe the dynamic behavior of the system accurately with small deformation assumption. Using the reduced-order model, the physical condition for the energy transfer was obtained by Hopf bifurcation analysis. From theoretical analysis, the critical values of parameters were obtained from the physical condition for the Hopf bifurcation.
The nonlinear behavior of the electrostatically actuated microbeam was analyzed via the pseudo-arclength continuation technique and a direct time integration. The frequency-amplitude responses and force-amplitude responses obtained by the pseudo-arclength continuation technique and a long-time-integration method showed good agreement. From the nonlinear electrodynamical resonance response, the mechanism of energy transfer between longitudinal vibration and transversal vibration was presented. The numerical simulations revealed that the bias voltage, the detuning frequency, and the excitation amplitude of in longitudinal direction played important roles in the capability of MEMS switches. Once the subcritical regime was activated, with a larger bias voltage or a smaller detuning frequency, shifted forced frequency responses were obtained, and the amplitude jump would become larger. The simulations gave a reliable analysis for reasonable design of structural properties, providing a great potential in improving precision of MEMS switches due to the subcritical bifurcation regime.