Mechanical Behaviors Research and the Structural Design of a Bipolar Electrostatic Actuation Microbeam Resonator

A class of bipolar electrostatically actuated micro-resonators is presented in this paper. Two parametric equations are proposed for changing the microbeam shape of the upper and lower sections. The mechanical properties of a micro-resonator can be enhanced by optimizing the two section parameters. The electrostatic force nonlinearity, neutral surface tension, and neutral surface bending are considered in the model. First, the theoretical results are verified with finite element results from COMSOL Multiphysics simulations. The influence of section variation on the electrostatic force, pull-in behaviors and safe working area of the micro-resonator are studied. Moreover, the impact of residual stress on pull-in voltage is discussed. The multi-scale method (MMS) is used to further study the vibration of the microbeam near equilibrium, and the relationship between the two section parameters of the microbeam under linear vibration was determined. The vibration amplitude and resonance frequency are investigated when the two section parameters satisfy the linear vibration. In order to research dynamic analysis under the case of large amplitude. The Simulink dynamics simulation was used to study the influence of section variation on the response frequency. It is found that electrostatic softening increases as the vibration amplitude increases. If the nonlinearity initially shows hardening behavior, the frequency response will shift from hardening to softening as the amplitude increases. The position of softening-hardening transition point decreases with the increase of residual stress. The relationship between DC voltage, section parameters, and softening-hardening transition points is presented. The accuracy of the results is verified using theoretical, numerical, and finite element methods.


Introduction
Nowadays, micro electromechanical system (MEMS) devices are attracting attention due to their small size, low weight, and low energy consumption. Among them, microbeams are a popular micro-component and are widely used in MEMS devices, e.g., as a microswitch [1,2], sensor [3,4], or resonator [5][6][7][8][9]. Therefore, designing and optimizing microbeams can improve the mechanical properties of MEMS devices. An electrostatically actuated microbeam resonator is presented in this paper.
Electrostatic actuation is one of the most common actuation methods in microstructures [10]. However, electrostatic actuation also has strong nonlinear characteristics [11]. The pull-in effect is the most critical unstable behavior [12]. When the actuation voltage exceeds a critical value, the microbeam will quickly come into contact with the plate, which can cause damage to the device in severe cases. Therefore, determining the unstable position of the microbeam, increasing the pull-in voltage, and the pull-in position of the resonator to improve its mechanical properties is a key issue in micro-resonator design. Abdel-Rahman et al. [13] presented a mathematical model for electrostatic actuation of microbeam resonators, which takes into account the neutral plane tension and nonlinear electrostatic force. It was found that the maximum dimensionless deflection reached 0.39 near pull-in, which is greater than the result calculated with the traditional method. Mobki, Harried et al. [14] presented two models describing unipolar plate of an actuated microbeam and an actuated bipolar plate microbeam; they found that the bipolar plate actuation model has a higher pull-in voltage. In addition, when the voltage ratio at both ends is 1, the increased gap distance ratio leads to decreased attraction voltage, and vice versa.
In addition to studying the pull-in behavior, vibration of the micro-resonator should also be analyzed. Nonlinearity can lead to multiple solutions, leading to jump instability. The relationship between the system's equivalent natural frequency and DC voltage depends on the strong nonlinearity and electrostatic softening. Moreover, the ideal working condition of the micro-resonator is linear vibration, thus the nonlinearity in the design should be eliminated to the greatest possible extent. Li et al. [15] studied the influence of physical parameters, such as microbeam viscoelasticity and DC voltage, on the system vibration. Finally, the parameters were optimized to determine the relationship between the physical parameters under linear vibration. Han et al. [16] used a spring-mass model to study the vibration of clamped-clamped microbeams and obtained a parameter design diagram relating the DC voltage and initial gap distance under linear conditions. The optimal design voltage under linear conditions was derived. In addition, nonlinearity can be utilized; Li et al. [17] used the modal coupling of a clamped-clamped microbeam to adjust the pull-in voltage and resonant frequency. Based on the coupled vibration behavior between the antisymmetric and symmetric modes, Hu et al. [18] presented a theoretical method for suppressing the midpoint displacement and reducing deformation. However, most early studies were based on rectangular microbeams and focused on optimizing their mechanical properties by adjusting their physical parameters.
As research has deepened, many scholars have found that the mechanical properties of a micro-resonator can be improved by adjusting the microbeam geometry when using the same raw materials. Joglekar and Trivedi [19,20] studied a class of clamped-clamped microbeams and proposed a parametric equation that describes smooth changes along the width of the microbeams. The system mechanical properties are improved by optimizing parameters in the relevant equation, and the results from the method were verified in several cases. On the basis of [19,20], Zhang et al. [21] discussed the influence of the optimized shape on the dynamic response of the microbeam. Kuang and Chen [22] optimized the mechanical properties of a micro-resonator by adjusting the microbeam thickness and gap distance. The required operating voltage range increased by a factor 6. Other studies focus on active control of micro-resonators; Alsaleem and Younis [23] introduced time-delay feedback into the dynamic control of electrostatically actuated micro-resonators. The influence of time delay displacement and velocity feedback on the system dynamic characteristics is discussed in terms of numerical and experimental results. Both kinds of control can improve the system dynamic behavior, enhance stability, and prevent pull-in under positive feedback gain. This is a very effective method for optimized resonator design.
So far, research on microbeam optimization primarily focuses on adjusting the microbeam size [24], using time delay feedback for active control [25,26], and optimizing the shape of the microbeam [27]. In this study, a micro-resonator was optimized by adjusting the microbeam shape, and the number of parametric equations was increased compared to models presented by predecessors. Equations are defined for the upper and lower sections of the microbeam such that the thickness changes uniformly along its length. Each equation has an independent parameter that can be used to adjust variations in the microbeam thickness. The resonator is optimized by choosing the appropriate values of two section parameters. The influence of section on the electrostatic force, small amplitude vibration, and dynamic behavior with large amplitude actuation in bipolar plates are investigated theoretically and with finite element simulations.
The structure of this paper is as follows. The mathematical model (partial differential equation) for electrostatic actuation of a microbeam resonator considering neutral stretching and bending is established in Section 2,. Two parametric equations are proposed for adjusting variations in the microbeam thickness. The Galerkin method is used to simplify the original vibration equation into an ordinary differential equation. In Section 3, the finite element method is used to simulate the electrostatic force on the system, and the influence of two section parameters on the electrostatic force, static pull-in behavior, and safe working area are discussed. In Section 4, the MMS is used to determine the response of the system under small amplitude vibrations. The relationship between the two section parameters is obtained by taking linear vibration as an optimization condition. Variation in the system's equivalent natural frequency and vibration amplitude are studied given this relationship. COMSOL Multiphysics was used to simulate and verify the natural frequency and electrostatic softening of the optimized micro-resonator. In Section 5, the influence of the section parameters on transitions between softening and hardening of the frequency response curve under large amplitude vibration is discussed. Moreover, the relationship between DC voltage, microbeam thickness, and the transition point is presented. Finally, a summary and conclusions are presented in the last section.

Governing Equation
A schematic of an electrically actuated microbeam is shown in Figure 1. The model consists of two fixed plates and a movable microbeam. Microbeams can be divided into six cases corresponding to various parametric equations. Figure 1a shows ordinary rectangular beams. Figure 1b,c shows a microbridge model, where the thickness does not change along beam length. Figure 1d,e shows thickened and thinned beams, respectively. both of which are symmetrical around the neutral plane. Figure 1f shows an irregular beam. A capacitor forms between the microbeam and the plate by applying a DC voltage V dc and an AC voltage V ac cos(Ωt) to the plate, thus generating an electrostatic force. and the number of parametric equations was increased compared to models presented by predecessors. Equations are defined for the upper and lower sections of the microbeam such that the thickness changes uniformly along its length. Each equation has an independent parameter that can be used to adjust variations in the microbeam thickness. The resonator is optimized by choosing the appropriate values of two section parameters. The influence of section on the electrostatic force, small amplitude vibration, and dynamic behavior with large amplitude actuation in bipolar plates are investigated theoretically and with finite element simulations. The structure of this paper is as follows. The mathematical model (partial differential equation) for electrostatic actuation of a microbeam resonator considering neutral stretching and bending is established in Section 2,. Two parametric equations are proposed for adjusting variations in the microbeam thickness. The Galerkin method is used to simplify the original vibration equation into an ordinary differential equation. In Section 3, the finite element method is used to simulate the electrostatic force on the system, and the influence of two section parameters on the electrostatic force, static pull-in behavior, and safe working area are discussed. In Section 4, the MMS is used to determine the response of the system under small amplitude vibrations. The relationship between the two section parameters is obtained by taking linear vibration as an optimization condition. Variation in the system's equivalent natural frequency and vibration amplitude are studied given this relationship. COMSOL Multiphysics was used to simulate and verify the natural frequency and electrostatic softening of the optimized micro-resonator. In Section 5, the influence of the section parameters on transitions between softening and hardening of the frequency response curve under large amplitude vibration is discussed. Moreover, the relationship between DC voltage, microbeam thickness, and the transition point is presented. Finally, a summary and conclusions are presented in the last section.

Governing Equation
A schematic of an electrically actuated microbeam is shown in Figure 1. The model consists of two fixed plates and a movable microbeam. Microbeams can be divided into six cases corresponding to various parametric equations. Figure 1a shows ordinary rectangular beams. Figure 1b,c shows a microbridge model, where the thickness does not change along beam length. Figure 1d,e shows thickened and thinned beams, respectively. both of which are symmetrical around the neutral plane. Figure 1f shows an irregular beam. A capacitor forms between the microbeam and the plate by applying a DC voltage dc V and an AC voltage cos( ) ac V t Ω to the plate, thus generating an electrostatic force.   Changes in the microbeam shape will produce an uneven electrostatic force between the upper and lower parts of the microbeam. The distance between the microbeam section at the fixed end and electrode plate is d, and the microbeam exhibits lateral motion under the action of electrostatic forces. A model whose thickness varies along microbeam length is considered here, where l and b are the length and the width of the microbeam, respectively. h is the thickness of microbeam at the clamped ends. The cross-sectional area and moment of inertia of the two clamped sides are A 0 = bh and I 0 = bh 3 /12, respectively. The microbeam thickness is determined by y 1 (x) = h 2 + λ 1 h sin πx l and y 2 (x) = − h 2 + λ 2 h sin πx l , where λ 1 and λ 2 are section parameters, and the microbeam section curvature changes with λ 1 and λ 2 . The cross-sectional area and moment of inertia are A(x) = A 0 (1 + (λ 1 − λ 2 ) sin(πx/l)) and I(x) = I 0 (1 + (λ 1 − λ 2 ) sin(πx/l)) 3 , respectively. In order to achieve better mechanical properties when both of microbeams consume same amount of material during their fabrication, it is assumed that the volume of the microbeam is constant. The relationship between h and λ 1 , λ 2 is found to be h = h 0 (π−2λ 1 +2λ 2 ) π by setting the microbeam volume does not change when the shape changes, where h 0 is the microbeam thickness when the beam is rectangular. The capacitance varies in microbeam different positions due to the changes in the gap between the microbeam and electrode plate. The capacitances of the upper and lower plates are where ε 0 is the dielectric constant in free space, ε r is the relative dielectric permittivity in the gap relative to free space, C 1 is the capacitance between the upper plate and the microbeam, and C 2 is the capacitance between the upper plate and the microbeam. The total energy stored by two capacitors is The electrostatic force is The equation of motion considers changes in the cross-sectional area and the electrostatic force on the bipolar plate, the equation is therefore rewritten as follows [28]: The boundary conditions are where c is the damping coefficient, and E and ρ are the elastic modulus and density, respectively. In Equation (4), N is the residual stress, the integral term is the tension and deformation of neutral plane.
For the convenience of calculation, consider following dimensionless: Substituting Equation (6) into (4) and (5) yields the dimensionless bending vibration equation: The corresponding boundary condition becomes where In the following figure, the " " symbol has been removed from the dimensionless parameters for convenience.
The fringing fields of electrostatic force are not taken into account in the above calculation. In order to observe the accuracy of electrostatic force, the electrostatic force in Equation (7) is compared with the electrostatic force considering edge effect. The results are shown in Figure 2. It can be found that the two results are consistent. Therefore, the infinite plate model is applied for convenience.
where c is the damping coefficient, and E and ρ are the elastic modulus and density, respectively.
In Equation (4), N is the residual stress, the integral term is the tension and deformation of neutral plane.
For the convenience of calculation, consider following dimensionless: Substituting Equation (6) into (4) and (5) yields the dimensionless bending vibration equation: The corresponding boundary condition becomes In the following figure, the "  " symbol has been removed from the dimensionless parameters for convenience. The fringing fields of electrostatic force are not taken into account in the above calculation. In order to observe the accuracy of electrostatic force, the electrostatic force in Equation (7) is compared with the electrostatic force considering edge effect. The results are shown in Figure 2. It can be found that the two results are consistent. Therefore, the infinite plate model is applied for convenience.

Galerkin Expansion
The deflection of the microbeam is defined as follows: The boundary conditions are u i (t) is the i-th modal coordinate amplitude, and φ i (x) is the i-th order undamped linear orthogonal mode function. According to results from the literature [29], a model with a single degree-of-freedom can sufficiently capture all key nonlinear aspects in the Galerkin approximation when studying the primary resonance. Assuming y(x, t) = u(t)φ(x), the existing conventional calculation is usually based on a Taylor expansion of the electrostatic force or by multiplying the denominator term of the electrostatic force [30]. The first method is a simple calculation, but the higher order displacement terms will be eliminated after Taylor expansion, which decreases the accuracy of the results. In order to accurately describe the electrostatic force while retaining the nonlinear characteristics of the electrostatic force as much as possible, a second calculation scheme is adopted in this paper. The denominator term is eliminated by multiplying (1 − y 1 (x) + h 2d − y(x, t)) 2 (1 + y 2 (x) + h 2d + y(x, t)) 2 at vibration equation. Since V ac is much smaller than V dc in the micro-resonator, the following calculation omits the ρ 2 term. Substituting Equation (9) into (7), multiplying by φ(x), and integrating from x = 0 to 1 yields the following equation: where, , v ac2 , and v ac3 are shown in Appendix A.
The static equilibrium equation can be written as follows: Without considering damping and AC disturbances, the following Hamiltonian system corresponding to the dimensionless dynamic equation can be obtained: The potential energy function is and the corresponding Hamiltonian is One should note here that the maximum lateral displacement of the microbeam is at the midpoint viz., The value of the modal function is φ(0.5) = 1.59, thus the range of u is u ∈ [ −1−λ 2 δ 1.59 , 1−λ 1 δ 1.59 ].

Static Analysis
First, a static analysis of the system is conducted. The safe region of the micro-resonator can be obtained from static analysis. Only by guaranteeing that the device can operate safely under the selected parameters can the research presented here be meaningful. Here, the accuracy of the electrostatic force obtained from theoretical analysis is verified. All possible microbeam cases are selected. The capacitance the between plates and microbeam was simulated using COMSOL Multiphysics, and the simulation model is shown in Figure 3. A voltage was applied at d on the clamped end of the microbeam while the microbeam was grounded. The entire model was placed in air. Assuming L = 400 µm, b = 45 µm, d = 3 µm, ρ = 2.33 × 10 3 kg/m 3 , E = 165 Gpa, and ε 0 = 8.85 × 10 −12 F/m, and the thickness of the clamped-clamped ends is h = 2 µm. The entire simulation was performed using the steady state solver. selected. The capacitance the between plates and microbeam was simulated using COMSOL Multiphysics, and the simulation model is shown in Figure 3. A voltage was applied at d on the clamped end of the microbeam while the microbeam was grounded. The entire model was placed in air. Assuming , and the thickness of the clamped-clamped ends is 2 μm h = . The entire simulation was performed using the steady state solver. The results from two methods are compared in Figure 4. This comparison shows that the two results are consistent. Although there is a slight deviation at larger voltage, the effect is insignificant. At this time, the voltage may have exceeded the pull-in voltage. Thus, high voltage cannot be used during operation. The specific pull-in voltage will be calculated later. Overall, the simulation results verify the correctness of theory.  The results from two methods are compared in Figure 4. This comparison shows that the two results are consistent. Although there is a slight deviation at larger voltage, the effect is insignificant. At this time, the voltage may have exceeded the pull-in voltage. Thus, high voltage cannot be used during operation. The specific pull-in voltage will be calculated later. Overall, the simulation results verify the correctness of theory. electrostatic force obtained from theoretical analysis is verified. All possible microbeam cases are selected. The capacitance the between plates and microbeam was simulated using COMSOL Multiphysics, and the simulation model is shown in Figure 3. A voltage was applied at d on the clamped end of the microbeam while the microbeam was grounded. The entire model was placed in air. Assuming , and the thickness of the clamped-clamped ends is 2 μm h = . The entire simulation was performed using the steady state solver. The results from two methods are compared in Figure 4. This comparison shows that the two results are consistent. Although there is a slight deviation at larger voltage, the effect is insignificant. At this time, the voltage may have exceeded the pull-in voltage. Thus, high voltage cannot be used during operation. The specific pull-in voltage will be calculated later. Overall, the simulation results verify the correctness of theory.  The theoretical results shown in Figure 4 are used to determine the influence of section variation on the electrostatic forces. The middle thickness of microbeam becomes thinner and the thickness of clamped-clamped ends becomes thicker when λ 1 is held constant and λ 2 decreases. The system potential can be reduced by decreasing λ 2 , as shown in Figure 5. This indicates that the potential energy is larger when the microbeam is thinner at the midpoint for a given voltage due to the greater electrostatic force. The theoretical results shown in Figure 4 are used to determine the influence of section variation on the electrostatic forces. The middle thickness of microbeam becomes thinner and the thickness of clamped-clamped ends becomes thicker when 1 λ is held constant and 2 λ decreases. The system potential can be reduced by decreasing 2 λ , as shown in Figure 5. This indicates that the potential energy is larger when the microbeam is thinner at the midpoint for a given voltage due to the greater electrostatic force. Next, the influence of the section parameters on static pull-in behaviors will be studied. In order to facilitate observation, the special case of a rectangular microbeam is added as a comparison, where the pull-in behavior is indicated with the black curve in Figure 6a . The microbeam is symmetric in the neutral plane in these two cases. Figure 6 shows that the static pull-in position gradually increases as the microbeam becomes thinner, eventually leading to secondary pull-in. Figure Figure 6a,b shows that the range of pull-in voltage will increase and secondary pull-in phenomenon will be more likely when the midpoint is thinner than the ends. In addition, static displacement jump will occur in the asymmetric microbeam. One can see from the case where 2 0.  Next, the influence of the section parameters on static pull-in behaviors will be studied. In order to facilitate observation, the special case of a rectangular microbeam is added as a comparison, where the pull-in behavior is indicated with the black curve in Figure 6a,b. The physical parameters of the microbeam are L = 400 µm, b = 45 µm, d = 3 µm, h 0 = 2 µm, E = 1.65 × 10 11 Pa, and ρ = 2.33 × 10 3 kg/m 3 . Figure 6a,b shows the case where the upper section parameters are λ 1 = 0.2 and λ 1 = −0.2. The lower section moves gradually from −0.3 to 0.3. At this time, the middle of the microbeam becomes thinner and the ends become thicker. The red curve corresponds to λ 1 = 0.2 and λ 2 = −0.2, while the green curve corresponds to λ 1 = −0.2 and λ 2 = 0.2. The microbeam is symmetric in the neutral plane in these two cases. Figure 6 shows that the static pull-in position gradually increases as the microbeam becomes thinner, eventually leading to secondary pull-in. Figure 6b shows the case where λ 2 = 0.2 and λ 2 = 0.3. Asymmetry of the microbeam will cause the static equilibrium point shift to the side with larger λ 1 and −λ 2 . Greater asymmetry leads to a larger migration distance. When λ 1 = −λ 2 , the microbeam is symmetric around the neutral plane and the equilibrium point position is symmetric around the x-axis. If λ 1 = −λ 2 > 0, the midpoint of the microbeam is thicker than at the ends, and vice versa if λ 1 = −λ 2 < 0. Comparing the three curves of red, black, and green in Figure 6a,b shows that the range of pull-in voltage will increase and secondary pull-in phenomenon will be more likely when the midpoint is thinner than the ends. In order to discuss this phenomenon more intuitively, Assume 82.5V . Similarly, in Figure 7b, is same as that of figure In addition, static displacement jump will occur in the asymmetric microbeam. One can see from the case where λ 2 = 0.3 in Figure 6b that the static displacement jumps from positive to negative when the DC voltage amplitude increases to V dc = 87 V. A comparison of Figure 6a,b shows that the static pull-in voltage mainly depends on the larger side of λ 1 and −λ 2 . As shown in Figure 6a, the static pull-in position of λ 2 from −0.3 to −0.2 varies greatly when λ 1 = 0.2, but the increasing trend in the static pull-in position decreases significantly when λ 2 ranges from −0.2 to 0.3. This arises because the static suction position at this time is mainly determined by λ 1 .
In order to discuss this phenomenon more intuitively, Assume V dc = 82.5V and draw the potential energy of λ 1 = −0.2 and λ 2 = 0.2 to observe the influence of another section variation on system safe area. When λ 1 = −0.2, the pull-in points of λ 2 = 0.19, λ 2 = 0.2, and λ 2 = 0.21 correspond to green, red, and blue points in Figure 7a, respectively. The lowest barrier energy is significantly lower when compared to the case when λ 2 = 0.   Residual stress may occur in the process of fabrication or operation of devices and the mechanical properties will be influenced. In order to accurately predict the safe working area of the device, the influence of residual stress on pull-in voltage and secondary pull-in will be discussed below. As shown in Figure 8, the pull-in voltage increases with the increase of residual stress. Figure  8a shows that when 1 =0.2 λ , the pull-in voltage increases with the increase of 2 λ and before Therefore, reducing λ 1 or increasing λ 2 in the design can improve the working range and stability of the resonator. On the other hand, static pull-in voltage primarily depends on the larger of λ 1 and −λ 2 . When one section of the microbeam is determined, the two fixed ends become thicker and more energy is lost at the ends when the other section is closer to the neutral plane, but the position of pull-in voltage does not change much. Therefore, taking λ 1 = −λ 2 can increase the pull-in voltage.
Residual stress may occur in the process of fabrication or operation of devices and the mechanical properties will be influenced. In order to accurately predict the safe working area of the device, the influence of residual stress on pull-in voltage and secondary pull-in will be discussed below. As shown in Figure 8, the pull-in voltage increases with the increase of residual stress. Figure 8a shows that when λ 1 = 0.2, the pull-in voltage increases with the increase of λ 2 and before λ 2 < −0.2, the pull-in voltage increases faster than after λ 2 > −0.2. This proves once again that the static pull-in voltage primarily depends on the larger of λ 1 and −λ 2 . Figure 8b shows the case of λ 2 = 0.2, It can be seen that the pull-in voltage decreases as λ 1 increases. Secondary pull-in phenonemon occurs in yellow and blue regions. Each λ 2 corresponds to three voltage values in this region. The dotted line indicates that the voltage on each line corresponds to the position in the static displacement picture. When λ 1 < −0.2, the secondary pull-in region is close to the lower plate and when λ 1 > −0.2, the secondary pull-in region is close to upper plate. In addition, the larger the yellow and blue areas, the larger the secondary pull-in range is. Therefore, the increase of residual stress will reduce the range of secondary pull-in area. λ > − , the secondary pull-in region is close to upper plate. In addition, the larger the yellow and blue areas, the larger the secondary pull-in range is. Therefore, the increase of residual stress will reduce the range of secondary pull-in area.

Dynamic Analysis
It is very important to conduct dynamic analysis of MEMS devices. The influence of section variation on the frequency response and bifurcation behavior can be observed through dynamic analysis, so as to optimize micro-resonator. MMS is used in the following to calculate small amplitude vibrations of the micro-resonator in the stable region. Here, we define First, D u can be obtained from Equation (12). Substituting D A u u u = + into Equation (11), ignoring the nonlinear terms above third order and nonlinear damping in Equation (11), and eliminating the term corresponding to the equilibrium position yields the small vibration equation at the equilibrium position: The meaning of each coefficient is shown in Appendix B. Because ac V is much smaller than dc V in the micro-resonator, we consider the cases In order to describe the primary resonance accurately, a tuning parameter σ is introduced and defined as The approximate solution to Equation (17) is expressed as follows:

Dynamic Analysis
It is very important to conduct dynamic analysis of MEMS devices. The influence of section variation on the frequency response and bifurcation behavior can be observed through dynamic analysis, so as to optimize micro-resonator. MMS is used in the following to calculate small amplitude vibrations of the micro-resonator in the stable region. Here, we define u = u D + u A , where u D is the response under DC voltage and u A is the response under AC voltage.
First, u D can be obtained from Equation (12). Substituting u = u D + u A into Equation (11), ignoring the nonlinear terms above third order and nonlinear damping in Equation (11), and eliminating the term corresponding to the equilibrium position yields the small vibration equation at the equilibrium position: ..
The meaning of each coefficient is shown in Appendix B. Because V ac is much smaller than V dc in the micro-resonator, we consider the cases V dc = O(1) and V ac = O(ε 3 ), where ε is a small dimensionless parameter. Equation (16) can be rewritten as ..
In order to describe the primary resonance accurately, a tuning parameter σ is introduced and defined as The approximate solution to Equation (17) is expressed as follows: where T n = ε n t, n = (0, 1, 2). Substituting Equations (18) and (19) into (17) and gathering like powers of ε yields: where D n = ∂ ∂T n , n = (0, 1, 2). The general solution to Equation (20) can be written as Substituting Equation (23) into (21) yields: where cc denotes the complex conjugate.
To eliminate the secular term, one needs which indicates that A is only a function of T 2 . Thus, Equation (24) becomes The solution u A2 can be written as follows: Substituting Equations (23) and (27) into Equation (22) yields the secular terms At this point, it is convenient to express A in polar form: Substituting Equation (29) into Equation (28) and separating the imaginary and real parts yields: where ϕ = σT 2 − β.
The steady-state response can be obtained by imposing the condition Da DT 2 = Dϕ DT 2 = 0. Finally, the frequency response equation can be derived as follows: where ξ = (10a q −a m ω 2 n )(a q −a m ω 2 n ) 24ω 3 n − 3(a c −a n ω 2 n ) 8ω n . The peak vibration amplitude and backbone curve are a max = f /(µω n ) and ω = ω n − κa max , respectively. The stability of the periodic solution can be determined using a method found in the literature [29,31].

Parameter Optimization
In the previous section, the stable region of the micro-resonator is optimized by static analysis, followed by conducting dynamic analysis. Previous research [29] has shown that the section parameter λ changes the softening and hardening behavior of the frequency response. It is noteworthy that there will be a moment of linear vibration in the transformation between softening and hardening, which is an ideal condition for micro-resonators. Therefore, the MEMS resonator is designed based on this condition. Assume d = 2 µm. First, one finds by calculation that the microbeam vibration is linear and rectangular when V dc = 23.95 V. The sectional parameters λ 1 and λ 2 under linear vibration are shown in Figure 9. respectively. The stability of the periodic solution can be determined using a method found in the literature [29,31].

Parameter Optimization
In the previous section, the stable region of the micro-resonator is optimized by static analysis, followed by conducting dynamic analysis. Previous research [29] has shown that the section parameter λ changes the softening and hardening behavior of the frequency response. It is noteworthy that there will be a moment of linear vibration in the transformation between softening and hardening, which is an ideal condition for micro-resonators. Therefore, the MEMS resonator is designed based on this condition. Assume 2 μm d = .
The coordinates of these special points listed in Table 1 are shown in Figure 9. The frequency response in nine cases was studied using MMS. In order to verify the above theoretical results, the long-time integral method (LTI) is used to calculate Equation (11) and the numerical solution is The coordinates of these special points listed in Table 1 are shown in Figure 9. The frequency response in nine cases was studied using MMS. In order to verify the above theoretical results, the long-time integral method (LTI) is used to calculate Equation (11) and the numerical solution is obtained. First, the frequency responses of D, H, and I are shown in Figure 10. Adjusting the AC voltage amplitude to a frequency response maximum amplitude of 0.1 shows that the vibration inside the curve exhibits hardening while the outer side exhibits softening. The point D corresponds to linear vibration in the special case of λ 1 = 0 and λ 2 = 0. obtained. First, the frequency responses of D, H, and I are shown in Figure 10. Adjusting the AC voltage amplitude to a frequency response maximum amplitude of 0.1 shows that the vibration inside the curve exhibits hardening while the outer side exhibits softening. The point D corresponds to linear vibration in the special case of 1 0 λ = and 2 0 λ = . Next, points A, B, C, and F in Figure 9 are analyzed to determine the influence of section variation on static pull-in so as to increase the working range of the micro-resonator, where point B and F are symmetric about y x = . Figure 11a shows that the maximum static suction voltage is  The equivalent frequency and amplitude of the microbeam in different conditions are also important for studying the micro-resonator. Figure 12 shows the frequency response curves in different conditions. It is found that the equivalent frequencies and amplitudes of various points in Figure 9 are different. In addition, the points that are symmetric about y x = in Figure 12 have the same equivalent frequency and different amplitudes. The vibration amplitude in region I is relatively large. Next, points A, B, C, and F in Figure 9 are analyzed to determine the influence of section variation on static pull-in so as to increase the working range of the micro-resonator, where point B and F are symmetric about y = x. Figure 11a shows that the maximum static suction voltage is obtained when λ 1 = −0.18 and λ 2 = 0.55. This also verifies the conclusion in Section 3, namely that the pull-in voltage is determined by the larger value out of λ 1 and −λ 2 . The static pull-in curve is symmetric around x = 0 when the point in Figure 7 is symmetric about y = x, where the static pull-in voltage is equal and pull-in position is the opposite, as shown in Figure 11b.  Next, points A, B, C, and F in Figure 9 are analyzed to determine the influence of section variation on static pull-in so as to increase the working range of the micro-resonator, where point B and F are symmetric about y x = . Figure 11a shows that the maximum static suction voltage is  The equivalent frequency and amplitude of the microbeam in different conditions are also important for studying the micro-resonator. Figure 12 shows the frequency response curves in different conditions. It is found that the equivalent frequencies and amplitudes of various points in Figure 9 are different. In addition, the points that are symmetric about y x = in Figure 12 have the same equivalent frequency and different amplitudes. The vibration amplitude in region I is relatively large. The equivalent frequency and amplitude of the microbeam in different conditions are also important for studying the micro-resonator. Figure 12 shows the frequency response curves in different conditions. It is found that the equivalent frequencies and amplitudes of various points in Figure 9 are different. In addition, the points that are symmetric about y = x in Figure 12    The variation in the resonance frequency and amplitude with the sectional parameters λ 1 and λ 2 are studied based on the results in Figure 9. Figure 13a shows that the resonance frequency is largest at point A, and the resonance frequency first decreases and then increases from point A to point B. From point B to point D, the resonance frequency decreases rapidly, reaching point λ 1 = 0, λ 2 = 0, and the resonance frequency is smallest. The trends of G, E, and F are the same as those of A, B, and C, and the corresponding equivalent frequency is equal at the point where y = x is symmetric. As shown in Figure 13b, the resonance amplitude increases from point A to point C. For the point corresponding to the y = x symmetry in Figure 9, the amplitude is larger when λ 1 is larger. The maximum amplitude appears between C and D.  The variation in the resonance frequency and amplitude with the sectional parameters 1 λ and 2 λ are studied based on the results in Figure 9. Figure 13a shows that the resonance frequency is largest at point A , and the resonance frequency first decreases and then increases from point A to point B . From point B to point D , the resonance frequency decreases rapidly, reaching point λ , and the resonance frequency is smallest. The trends of G, E , and F are the same as those of A , B , and C , and the corresponding equivalent frequency is equal at the point where

Dynamic Analysis with Large Amplitude
Since traditional MMS is only suitable for small vibrations, the traditional MMS will produce incorrect results when vibration amplitude exceeds a certain value [32]. In order to study the mechanical properties of the system when the vibration exceeds a certain amplitude, a new multi-scale method (NMMS) is introduced by combining the homotropy idea with the multi-scale method [33]. Equation (11) is solved using a Simulink dynamics simulation in order to verify the accuracy of the results.
It is assumed that the final vibration frequency is equal to the excitation frequency. A scaling parameters ε can be used to convert Equation (11) into the following form ..
The approximate solution to Equation (33) is expressed as follows: Substituting Equation (34) into (33) and equating coefficients of like powers of ε yields the following equations: where: D n = ∂ ∂T n , n = (0, 1, 2). The general solution to Equation (35) can be written in the following form: where A(T 1 ) is vibration amplitude and β(T 1 ) is the vibration phase. The secular term derived from MMS is expressed as: where ϕ = ωT 0 + β(T 1 ). Substituting Equation (37) into Equation (38) and separating the imaginary and real parts yields: The steady-state response can be obtained by applying the condition Da/DT 1 = Dϕ/DT = 0. The ultimate frequency response can be deduced as follows: The micro-resonator vibration amplitude increases as the AC voltage amplitude increases. The physical parameters of the microbeam are L = 400 µm, b = 45 µm, d = 3 µm, h 0 = 2 µm, V dc = 60 V, E = 1.65 × 10 11 Pa, and ρ = 2.33 × 10 3 kg/m 3 . Select the upper section parameters λ 1 = −0.2 and λ 2 = 0.2. The results from MMS, NMMS, and the Simulink dynamics simulation are compared in Figure 14, and the results from the three methods are found to be consistent when V ac = 0.05 V. However, the MMS results begins to exhibit errors when V ac = 0.1 V as the vibration amplitude increases, while the NMMS results are still consistent with the simulation results. The system vibration appears to transition from hardening to softening when V ac = 0.2V. The new behavior can be obtained with NMMS. It is worth noting that the frequency response curve is composed of high and low energy branches when the vibration amplitude is small. As the amplitude increases, the two branches gradually converge. Since the main purpose of this chapter is to discuss the nonlinear softening-hardening transition behavior, it will not be discussed in further detail here. A detailed discussion can be found in literature [28].  Two kinds of motion can occur at the transition point: dynamic pull-in and jump to a higher stable branch. The energy output is higher in the second case. Changes in the softening-hardening transition point can be observed as the vibration amplitude increases from V ac = 0.2 V to V ac = 0.25 V, V ac = 0.3 V, and V ac = 0.35 V, and the corresponding frequency response curve can be drawn. One can see from Figure 15 that the displacements corresponding to nonlinear softening-hardening transition points in the four cases are almost the same. It can be inferred that the position of the transition point does not change as the vibration amplitude increases or the change is negligible at the transition point.
According to previous studies, the location of the transition point is related to the microbeam thickness and DC voltage. The location of the transition point can be obtained by calculating du/dω = 0, yielding the relationship between the thickness, DC voltage, and transition point. The transition point can only be obtained numerically because the resulting equation is more complicated. Curve fitting is used here to determine the transition point. There are two transition points in each frequency response curve. Figure 16 shows that the amplitudes corresponding to the two transition points are nearly equal. Therefore, only one curve can be used for the two positions. As can be seen from Figure 14, the position of softening-hardening transition point decreases with the increase of residual stress. Figure 16a shows the relationship between the thickness and the transition point position, and Figure 16b shows the relationship between the DC voltage and the transition point position. The straight line shows the solution drawn with the fitted equation, and the circle shows the numerical solution obtained by the Newton iteration. The green and blue circles show the left and right transition points, respectively. A comprehensive analysis of the two curves yields the following equation:  Curve fitting is used here to determine the transition point. There are two transition points in each frequency response curve. Figure 16 shows that the amplitudes corresponding to the two transition points are nearly equal. Therefore, only one curve can be used for the two positions. As can be seen from Figure 14, the position of softening-hardening transition point decreases with the increase of residual stress. Figure 16a shows the relationship between the thickness and the transition point position, and Figure 16b shows the relationship between the DC voltage and the transition point position. The straight line shows the solution drawn with the fitted equation, and the circle shows  Curve fitting is used here to determine the transition point. There are two transition points in each frequency response curve. Figure 16 shows that the amplitudes corresponding to the two transition points are nearly equal. Therefore, only one curve can be used for the two positions. As can be seen from Figure 14, the position of softening-hardening transition point decreases with the increase of residual stress. Figure 16a shows the relationship between the thickness and the transition point position, and Figure 16b shows the relationship between the DC voltage and the transition point position. The straight line shows the solution drawn with the fitted equation, and the circle shows Next, the relationship between the DC voltage, microbeam thickness, and transition point position is verified. Table 2 shows the case where the voltages are V dc = 30 V, V dc = 40 V, V dc = 50 V, and V dc = 60 V. Each case considers four different thicknesses to account for the breadth of the data. u e is the transition point position obtained from empirical equation, u n l and u n r are the positions of left and right transition points obtained by setting du/dω = 0, respectively, u n l is the left transition point, and u n r is the right transition point. The "/" symbol indicates no transition point, i.e., there is no softening-hardening transition. Table 2 shows that the error between the transition point position and the actual transition point obtained from Equation (41) is small, which verifies the accuracy of the equation.

Equivalent Natural Frequency Simulation
The microbeam resonator was explored and optimized using theoretical and numerical analysis. The finite element method will be used to simulate the optimization results using COMSOL Multiphysics software.
In COMSOL, the "solid mechanics" interface and "electrostatic" interface are combined with a dynamic grid function. First, the steady state of the model and the eigenfrequency in the steady state are studied. Here, finite element simulations are performed for the aforementioned seven cases, and the calculated simulation results are shown in Figure 17. One can see that the theoretical and simulation results are consistent, and the pull-in voltage at A, B, C, D, and G in Figure 17 are consistent with those in Figure 11.
The seven cases correspond to linear vibration when V dc = 23.95 V. The equivalent theoretical frequency results are compared with the finite element equivalent frequency results when V dc = 23.95 V in Table 3. It is found that the maximum error is only 5.49%, which shows that the theoretical and numerical methods are consistent. are studied. Here, finite element simulations are performed for the aforementioned seven cases, and the calculated simulation results are shown in Figure 17. One can see that the theoretical and simulation results are consistent, and the pull-in voltage at A, B, C, D, and G in Figure 17 are consistent with those in Figure 11.  Table 3. It is found that the maximum error is only 5.49%, which shows that the theoretical and numerical methods are consistent.

Equivalent Natural Frequency Simulation
The simulated equivalent natural frequency contains some error, although the error between the finite element solution and the theoretical solution is relatively small. More simulated frequency response curves are presented in this section. The physical parameters of the microbeam are L = 400 µm, b = 45 µm, d = 3 µm, h 0 = 2 µm, V dc = 60 V, E = 1.65 × 10 11 Pa, and ρ = 2.33 × 10 3 kg/m 3 . The simulated sectional parameters were λ 1 = −0.2 and λ 2 = 0.2. Figure 16 shows a comparison between the simulation and theoretical results.
One can see that the amplitude response is linear when the amplitude is small, while hardening appears at higher amplitude. The red asterisk shows the results obtained from the frequency response analysis in COMSOL. One can see that the simulated amplitude is consistent with the theoretical amplitude when the AC voltage varies. When V ac = 0.02V, the simulated amplitude is 0.0409 and the theoretical amplitude is 0.0406, yielding a relative error of −0.73%. When V ac = 0.1 V, the simulated amplitude is 0.204 and the theoretical amplitude is 0.203, yielding a relative error of 0.49%. However, because the frequency response analysis solver is a linear solver, the results do not exhibit hardening or softening. Each frequency is selected separately in order to further simulate the nonlinearity. After determining the frequency at a certain point, the corresponding amplitude is obtained with the transient solver. This step was repeated with different frequencies in order to determine the corresponding amplitude. This result is shown in blue circles in Figure 18. One can see that the finite element solution exhibits hardening. Moreover, the simulated natural frequency in the frequency domain is the same as the natural frequency in the time domain. The two simulation results are very consistent when the amplitude is small, but the maximum vibration amplitude from the time domain is lower than the other two results as the amplitude increases. This phenomenon may arise because the time interval corresponding to the highest point position is too short in the time domain analysis, thus the displacement does not point near the peak. The simulated equivalent natural frequency contains some error, although the error between the finite element solution and the theoretical solution is relatively small. More simulated frequency response curves are presented in this section. The physical parameters of the microbeam are  Figure 16 shows a comparison between the simulation and theoretical results. One can see that the amplitude response is linear when the amplitude is small, while hardening appears at higher amplitude. The red asterisk shows the results obtained from the frequency response analysis in COMSOL. One can see that the simulated amplitude is consistent with the theoretical amplitude when the AC voltage varies. When 0.02V ac V = , the simulated amplitude is 0.0409 and the theoretical amplitude is 0.0406, yielding a relative error of −0.73%. When 0.1 V ac V = , the simulated amplitude is 0.204 and the theoretical amplitude is 0.203, yielding a relative error of 0.49%. However, because the frequency response analysis solver is a linear solver, the results do not exhibit hardening or softening. Each frequency is selected separately in order to further simulate the nonlinearity. After determining the frequency at a certain point, the corresponding amplitude is obtained with the transient solver. This step was repeated with different frequencies in order to determine the corresponding amplitude. This result is shown in blue circles in Figure 18. One can see that the finite element solution exhibits hardening. Moreover, the simulated natural frequency in the frequency domain is the same as the natural frequency in the time domain. The two simulation results are very consistent when the amplitude is small, but the maximum vibration amplitude from the time domain is lower than the other two results as the amplitude increases. This phenomenon may arise because the time interval corresponding to the highest point position is too short in the time domain analysis, thus the displacement does not point near the peak.

Conclusion
Static and dynamic analyses of electrostatically-driven microbeam resonators are presented in this paper. Two parametric equations describing control over the upper and lower sections of microbeam were introduced. The mechanical properties of the micro-resonator are strengthened by optimizing the parameters in these equations. The neutral plane tension, neutral plane bending, and electrostatic nonlinearity are considered in the mathematical model. First, the model is reduced and simplified by Galerkin discretization. Then the relationship between two section parameters under

Conclusions
Static and dynamic analyses of electrostatically-driven microbeam resonators are presented in this paper. Two parametric equations describing control over the upper and lower sections of microbeam were introduced. The mechanical properties of the micro-resonator are strengthened by optimizing the parameters in these equations. The neutral plane tension, neutral plane bending, and electrostatic nonlinearity are considered in the mathematical model. First, the model is reduced and simplified by Galerkin discretization. Then the relationship between two section parameters under linear vibration is obtained by a Newton iteration. Finally, the vibration at small and large amplitudes were obtained using MMS and NMMS. The results were verified with the Runge-Kutta method, Simulink dynamic simulation and finite element method.
The key conclusions are as follows.
1. The potential energy and electrostatic force are greater when the microbeam intermediate thickness is thinner for a given applied voltage. Decreasing the thickness at the midpoint of the microbeam will increase the pull-in voltage and increase the likelihood of secondary pull-in. The pull-in position is biased toward the larger of λ 1 and −λ 2 when microbeam is asymmetric about the neutral plane, and the pull-in voltage is primarily determined by the larger of λ 1 and −λ 2 .

2.
When the micro-resonator vibrates linearly, the relationship between λ 1 and λ 2 is symmetric at y = x. The eigenfrequency and amplitude of the microbeam are studied when λ 1 and λ 2 correspond to linear vibration. Regarding the y = x symmetric point in Figure 7, the microbeams have the same eigenfrequency and different amplitudes. The amplitude near the AC voltage source is larger.

3.
The MMS should not be used as the amplitude increases. The frequency response obtained with MMS will exhibit weaker nonlinear softening. If the initial vibration exhibits hardening, the transition between hardening and softening will occur at a certain time as the vibration amplitude increases. The position of the transition point does not change as the amplitude continues to increase. At the transition point, the micro-resonator can jump to a higher stable branch and output more energy.
It is found that the static and dynamic behaviors of micro-resonator can be changed by adjusting the microbeam shape. Determining the appropriate geometry can improve the mechanical properties for given raw materials.

Conflicts of Interest:
The authors declare no conflict of interest.