Dynamic Characteristics of Electrostatically Actuated Shape Optimized Variable Geometry Microbeam

We mainly analyze the dynamic characteristics of electrostatically actuated shape optimized variable geometry microbeam. A nonlinear dynamic model considering midplane stretching, electrostatic force, and electrical field fringing effects is developed. Firstly, we study the static responses of the optimized microbeams under DC polarization voltage. The generalized differential quadrature method (GDQM) is used. Secondly, the dynamic responses of the shape optimized microbeams driven by DC and AC voltages are investigated using GDQM in conjunction with Levenberg-Marquardt optimization method. The results show that the more gradual change inwidth, the larger the resonant frequency and themaximumamplitude at resonance.Thenwe further discuss in detail how do the maximum width, midsection width, and curvature of the width function affect the frequency response of the microbeams. We find that the amplitude and resonant frequency of the dynamic response are not monotonically increasing as the curvature of the width function increases and there exists a critical curvature. This analysis will be helpful in the optimal design of MEMS actuators. Finally, for more consideration, different residual stress, squeeze-film damping, and fringing effect models are introduced into the governing equation of motion and we compare the corresponding dynamic response.


Introduction
Electrostatically actuated microbeams have been widely used in microelectromechanical systems (MEMS) devices such as capacitive MEMS switches [1] and gyroscope [2].An electrostatically actuated microbeam consists of a movable electrode (microbeam) and a fixed electrode (substrate), with dielectric medium filling the gap between them (Figure 1).This structure is simple and well compatible with existing microfabrication techniques.DC and AC voltages are combined to drive the microbeam.The AC voltage value is much smaller than that of the DC voltage.The DC voltage deflects the microbeam to a new equilibrium position, while the AC voltage creates oscillatory motion around this equilibrium position.
However, electrostatic actuation has an intrinsic instability situation, known as the pull-in phenomenon [3,4].When the applied voltage exceeds a critical value, which is called the pull-in voltage, the mechanical restoring force can no longer resist the electrostatic force, thereby causing the collapse of the microbeam.In most cases, pull-in instability greatly limits the stable operation range of the microbeam.However, a large stable range can be extremely useful in a wide variety of tuning applications.Many scientists therefore have done much research to expand the stable range.A comprehensive comparison of various methods of stable range extension has been done by Zhang et al. [5].
Recently, researchers have paid their attention to variable geometry microbeams [6].Since the microbeam shape can influence both the structural stiffness and the electrostatic force, it would be a very effective method to improve the performance of the microdevices by optimizing the microbeam shape.More and more researchers have been devoted to the optimization work, and the main optimization methods are divided into two categories: topology optimization and continuous geometry optimization.
Sigmund [7] studied an electrothermally driven microactuator with topology optimization method.They discussed the differences between modeling and optimizing the actuators using linear and nonlinear finite element analyses.Raulli and Maute [8] analyzed the design of electrostatically actuated MEMS by topology optimization.They revealed the advantages of varying simultaneously the interface topology and the layout of the electrode versus conventional approaches of optimizing only the structural layout.Lemaire et al. [9] described how to utilize the topology optimization in the multiphysics field of MEMS to maximize pull-in voltage of an electromechanical device.Although these studies have indeed improved the performance of the microdevice, the resulting optimized structure using topology optimization contains some closed holes, which is hard to realize with the existing fabrication technology.Therefore, several researchers turn to explore continuous geometry optimization.Abdalla et al. [10] proposed a passive method to maximize the pull-in voltage for an electrostatically actuated microbeam by changing the microbeam shape in width and thickness.Their results illustrated that width optimization seems to be more beneficial than thickness optimization.However, they ignored the fringing effects at the electrode edges and the effects of geometric nonlinearity.Najar et al. [11,12] investigated the static and dynamic behavior of a variable section microactuator.They first demonstrated the convergence of the differential quadrature method (DQM) in conjunction with the finite difference method (FDM) and successfully used this method to analyze the effects of variations in the gap size and microbeam width and thickness on the frequency response.But they did not consider the fringing effect either.Joglekar and Pawaskar [13] proposed a versatile parametric width function which smoothly varies the width of the fixed-fixed microbeam along its length.They optimized the parameters of the proposed width function and enumerated several cases to demonstrate their methodology.The goal of their research is to enhance the static and dynamic pull-in ranges of electrostatically actuated microbeams, but they did not discuss the impact of the resulting optimized shape on the dynamic response of the microbeam.Herrera-May et al. [14][15][16] not only studied the resonant characteristic of the single layer variable cross-section microbeam but also researched the bending resonant frequency of multilayered microresonators with variable cross-section.So far, the research on the dynamic characteristic of the optimized microbeam is fewer.In this paper, we carry out the dynamic analysis of the optimized microbeam and discuss the impact of adjusting the geometric parameters on the frequency response of the shaped microbeam.
This paper is organized as follows: in Section 2, we introduce the mathematical model, accounting for the system nonlinearities due to midplane stretching, electrostatic force, and fringing effects.Then in Section 3, the static pull-in parameters (voltage and travel range) and dynamic characteristics of optimized microbeams are extracted.The results of our calculation are compared with those reported in the literature.The influences of varying the axial load, fringing effect formulas, and microbeam shape on the dynamic response of the optimized microbeam are shown in Section 4. Finally, Section 5 is the conclusion of the paper.

Mathematical Model
In this paper, we consider a variable cross-section fixed-fixed microbeam (Figure 1) driven by a combination of DC and AC voltages.The microbeam is suspended above a fixed electrode with an initial gap  0 .Both the movable and fixed electrodes are good conductors of electricity and the permittivity of free space between them is .As combined driving voltages are applied between the movable and the fixed electrodes, a position-dependent electrostatic force is distributed and deforms the movable electrode towards the fixed electrode.
In order to consider the influence of the fringing effects, an extra term is introduced into the governing equation [12].The partial differential integral equation for the transverse deflection  is written as (1)

Shock and Vibration 3
The boundary conditions are given by  (0, ) = 0,   (0, ) = 0,  (, ) = 0,   (, ) = 0, where  is the material density, () and () represent the cross-section area and the area moment of inertia,  is the effective Young's modulus, and   and V() are, respectively, the DC and AC driving voltages.Equation (1) considers residual stress through the parameter .The geometric nonlinearity due to midplane stretching is also taken into consideration by means of integral differential term.The first term on the right side of (1) represents the electrostatic force considering the fringing effects.
The following dimensionless variables and parameters are introduced: where  0 and ℎ 0 represent, respectively, the width of the microbeam at the clamped edge and the microbeam thickness.b() is the variable width function along the axial direction.Substituting (3) into (1) and (2), the dimensionless form of the governing equation is rewritten as and the associated boundary conditions become The parameters in (4) are defined as follows:

Solution Method
Considering the nonlinearities present in (4), the GDQM [17] One essential issue pertaining to this method is how to efficiently and accurately choose the weighting coefficients    .In general, a poor choice of such weights leads to a numerical ill position of the problem.Based on the theory of Lagrange interpolation, the GDQM uses a Lagrangian interpolation polynomial in determining the weighting coefficients: V  = 1, 2, . . .,  = .
The weighting coefficients of higher order derivatives can be determined by a recurrence relationship: where  is the total number of discrete grid points.In order to have more accurate solutions, the Chebyshev-Gauss-Lobatto distribution is adopted: The integral terms in the governing equation are computed by the Newton-Cotes formula and the cotes coefficients are computed by

Static Solution.
Neglecting all the time-dependent terms in (4), the governing equation of static behavior of the microbeam is given by with In the following, three cases of prismatic microbeams are considered.We calculate the pull-in parameters and the comparison with the previous reference in literatures is listed in Table 1.The results of present work shown in Table 1 closely agree with those reported in the literature, which verifies the accuracy of the method.
Figure 2 shows the maximum deflection-voltage relationship of the third case in Table 1, from zero up to the pull-in voltage.The instability occurs when the slope of the deflection-voltage curve becomes infinite.Our results are in good agreement with those in [19], where they solved a reduced order model [20] using Galerkin method.
In order to enhance the static pull-in travel range, Joglekar and Pawaskar [13] proposed a novel versatile parametric and continuous width function.They use energy-based techniques to extract static pull-in parameters and use Nelder-Mead method to optimize the parameters of the proposed width function.Nine groups of optimized parameters were put forward and it was pointed out that these nine cases are applicable to any microbeam geometry having the ratio of  0 / = 0.15.This paper uses the above method to compute the static pull-in parameters and compares the results with those reported in the literature.The comparison is listed in Table 2. Above the oblique line are the results of our calculation, while under the oblique line are the results reported in the literature.
Compared with the prismatic microbeam, the pull-in travel ranges are greatly improved after the shape optimization (see Table 2).Among all the nine cases, the third one improved the most (up to 19.66%).Some of the cases have relative larger differences between our results and the reference.This is caused by the approximate method we adopted when dealing with the optimized microbeam with sharp points (the profile is shown in next section).Because the GDQM is based on the theory of Lagrange interpolation, it has some drawbacks when handling the microbeam with derivative discontinuous width function profile.To improve the accuracy, we approximate the width function by high order polynomial, and this approximation will introduce several errors.Nevertheless, the results are acceptable.
The analysis above only illustrates the maximum deflection of the microbeam with voltages.It cannot intuitively reflect the deflection of the whole microbeam.In order to understand the overall deformation of the microbeam when the driven voltage is close to the pull-in voltage, we also extract the displacement of the other discrete points besides the midpoint of the microbeam.Figure 3 displays the overall deformation of the nine shape optimized microbeams.
3.2.Dynamic Solution.Some shape optimization results are given in the literature [13] and these optimized microbeams indeed greatly increase the travel range.But the literature did not discuss such optimization shapes would have what kind of influence on the dynamic response of the microbeam.Based on their research, we further analyze the dynamic behavior of the shape optimization microbeams.For convenience, we rewrite the dimensionless governing equation (4) as the following form.The prime denotes the derivatives with respect to the space coordinate  and the overdot denotes the derivatives with respect to time : Newton-Cotes formula and GDQM are applied, respectively, to discretize the integral and derivative terms in (14).Then the time of one period is discretized into  points.In order to solve the obtained ordinary differential equations, finite differential method is used to compute the velocity and the acceleration of the microbeam and subsequently make the results satisfy (15); finally utilizing the Levenberg-Marquardt optimization method, the deflection  with respect to  and  can be obtained: where   = 2/Ω( − 1), Ω is the excitation frequency  , = (  ,   ), and  V , =  V (  ,   ) = ẇ (  ,   ).Using this method, we calculate a typical example of microbeam proposed in [19] and obtain its dynamic behavior.Figure 4 indicates that our estimates of the dynamic response in both hardening and softening domains closely agree with those reported in the literature, where Nayfeh and coworkers got the similar results by a combination of two-point boundary-value problem and shooting method [21].This comparison verifies the accuracy of the method we adopted.
As described in literature [22], Younis and Nayfeh found that DC voltage affects the qualitative and quantitative nature of the effective nonlinearity coefficient of the system, which changes the dynamic response of the microbeam from a hardening to a softening response.This is because, with increasing of the DC voltage, the electric nonlinearity, which tends to lower the response frequency, drastically dominates the influence of the geometric nonlinearity, which tends to increase the resonant frequency.They concluded that most of the models used in the literature neglect the effect of the electric nonlinearity and particularly the quadratic one.Usually, these models assume the nonlinearity of the system to be solely cubic and positive, which predicts a hardening behavior rather than the correct softening behavior.
In the following sections, we first assigned  dc = 23 V and got the softening curves and then changed the DC voltage to  dc = 10 V, getting the hardening curves.
Next, we adopt this method to analyze the dynamic response of the shape optimization microbeams proposed in [13].The parametric width function (first put forward by Joglekar and Pawaskar in [23]) for fixed-fixed microbeam is shown as follows: and the acceptable ranges of the four parameters are where  min and  max are constraints on the minimum and maximum allowable width and  0 is the width of the referential prismatic microbeam.Table 3 shows the optimized parameters of dynamic mode of actuation of fixed-fixed microbeam.The main reason for choosing dynamic mode parameters instead of static ones is that our analysis target is dynamic response.In addition, most of the optimization shapes of the static and dynamic modes closely resemble each other.For more accuracy, the optimized results of dynamic mode (Table 3) are chosen as the basis of calculation.Table 3 is divided into three groups: D1-D3, D4-D6, and D7-D9.Each group has the same minimum constraints  min at the midpoint of the microbeam.All the nine cases are computed and the calculated frequency response curves are shown in Figures 5-7.
The green curve in Figure 5 stands for the frequency response of the reference prismatic microbeam, and the other three are those of the optimized microbeams: from left to right are in turn Case 1 to Case 3. The rectangular box displays the schematic diagram of the optimized width shapes of the microbeam.The shape of the referential prismatic microbeam is also drawn on each plot in order to facilitate a visual comparison between the referential and the optimized shape.All these three optimized microbeams have the same minimum width value at the midpoint.Compared with the prismatic microbeam, the amplitude of the frequency response of Case 3 improves the most; that is to say, the geometry shape of Case 3 is better than the other two.By comparing the three different shapes, we can conclude that the more gradual change in width, the larger the maximum amplitude at resonance.And this conclusion is also applicable to the following two groups.As can be seen from Figure 5, the microbeam shape of Case 1 contains obvious stress concentration phenomenon.While in practical application, such a shape of the microbeam always avoids being manufactured.
A constant volume constraint is taken into account during the optimization procedure, which helps in comparing the prismatic and optimized microbeams for their performance.All the microbeams consume same amount of material during their fabrication.The different shapes of the microbeam can be seen as a new distribution of the material.It can be observed from Figure 6 that the midsection of Case 6 is longer than that of Case 5; meanwhile, the resonance frequency and amplitude of Case 6 are also larger.Through the comparison, it can be deduced that the longer the constant width portion near the central, the larger the amplitude of the frequency response.
The optimization parameters of Case 8 and Case 9 are identical; thus, the frequency responses overlap together.We can see from Figures 5-7 that the amplitude and resonance frequency of the optimized microbeams are all increased compared to the reference prismatic microbeam.Among the nine optimal shapes, the third one is the best, which agrees with the results of the static solution.The shape of Case 1, Case 4, and Case 7 all contains stress concentration parts, and these structures only make sense in mathematical analysis and they will never be put into practical application.

Numerical Analysis
There are four variable parameters in the width function (16).
Different parameters have different effects on the microbeam shape, and the changes in the shape will further affect the dynamic behavior of the microbeam.Joglekar et al. put forward a set of optimized parameters, but they did not discuss the influence of each parameter on the dynamic response of the microbeam.In this section, we emphasize on studying the influence of different shape parameter changes of the microbeam on the dynamic response.Then, for comprehensive consideration, we also consider the impact of different axial loads, squeeze-film damping, and different fringing effect models on the frequency response.Without loss of generality, a typical example of fixed-fixed microbeam proposed in [13] has been selected.The physical parameters, material properties, and constraints are listed in Table 4. Figure 8 depicts the optimized width shape of this microbeam as well as the referential prismatic microbeam.

Geometry Effects.
This subsection discusses the effect of varying the geometric shapes on the microbeam frequency response.The maximum width (), midsection width (, ), and curvature () are varied.We study the dynamic behavior of the microbeam under  dc = 23 V and  ac = 0.5 V.The results are shown, respectively, in Figures 9-13.
The acceptable ranges of  are  min / 0 ≤  ≤  max / 0 .As for this microbeam, 0.4 ≤  ≤ 2. Figure 9 shows the frequency responses with respect to different .We choose  four values of , and the difference between each of them is 0.5.In the bottom left corner of Figure 9, a schematic diagram of the optimized microbeams with different  is displayed.While at the top right corner is the partial enlargement of the response curves.It can be seen that with the increase of , the response amplitude increases.However, the discrepancy between the curves is very small.This is because augmenting the maximum width () not only renders the microbeam stiffer but also increases the electrostatic force due to the increase of the overlap area of the electrode plates.It can be found that the discrepancy of the response tends to be smaller with the increase of .

Shock and Vibration
Parameters  and  control the midsection width of the microbeam.We can see from Figures 10 and 11 that increasing the value of () decreases (increases) the resonance frequency and maximum amplitude at resonance.From the perspective of physical analysis, increasing  (decreasing ) not only expands the width of the microbeam (enhancing the stiffness and mechanical restoring force), but also augments the overlap area of the capacitor (enhancing the electrostatic force).With this kind of shape change, the narrower the microbeam, the smaller the equivalent stiffness and mass of the microbeam, but the descent velocity of the mass is larger than that of the stiffness; thus, the natural frequency becomes larger as the microbeam width narrows down.When expanding the frequency sweep range, we also find superharmonic resonance phenomenon (Figure 10).The superharmonic resonance shows almost linear behavior and this agrees with the conclusions of [24].
Changes of  do not change the width of the microbeam at the midpoint.Augmenting  increases the curvature of the microbeam's width curve.When  is small, augmenting  increases the resonance frequency as well as the maximum amplitude at resonance (Figure 12).However, when  continues to increase (larger than 1.4), the trend becomes just the opposite, the response begins to shift to the lower frequency, and the maximum amplitude also begins to decrease, which is shown in Figure 13.This is because when  is smaller than 1.4, increasing  renders the descent velocity of the equivalent stiffness of the microbeam smaller than that of mass, just like the change regularity of parameter .When  is larger than the critical value 1.4, as the increasing of , the rate of decline of equivalent stiffness of the microbeam becomes faster than the microbeam mass, thus obtaining the opposite change trend.This kind of shape change is very interesting and it is beneficial to the optimal design of MEMS resonators because of the wide adjustable range.

Residual Stress Effects.
During the process of micromanufacturing, it is inevitable to produce residual stress along the axial direction of the microbeam.The stress can be tensile or compressive, which will accordingly increase or decrease the microsystem structural stiffness.To ensure accuracy, we take into account this effect through the parameter .
Midplane stretching is by far the most important source of geometric nonlinearity in MEMS and it increases the structural stiffness of the microbeam in a nonlinear way that resembles a cubic effect [21].In contrast, electrostatic force softens the modal stiffness of the structure.In addition, the softening effect gets stronger with the increase of the DC voltage.So when the driving voltage is large, the influence of midplane stretching is weaker than electrostatic force; thus, all the frequency response curves in Figures 5,6,7,9,10,11,12,and 13 show softening behavior.For more comprehensive consideration, in this subsection, we choose a smaller DC voltage  dc = 10 V; meanwhile, AC voltage remains constant.In this case, the softening effect caused by the electrostatic force becomes weaker, so the system has a hardening effective nonlinearity [24].The microbeam shape is still the optimization shape shown in Figure 8.
Figure 14 shows the frequency response of the midpoint of the fixed-fixed optimized microbeam for five different cases of axial load.The  in Figure 14 is a nondimensional quantity, positive value represents the axial tensile stress, and negative value represents the axial compressive stress.It can be seen from Figure 14 that augmenting the tensile stress increases the resonance frequency and reduces the maximum amplitude at resonance, while augmenting the compressive stress reduces the resonance frequency and increases the maximum amplitude at resonance.Here also find superharmonic resonance phenomenon as outlined in Figure 14.

Squeeze-Film Damping Effects.
In MEMS, due to the scale effect of the microstructure, air damping has obvious influence on the dynamic behavior of the microdevice.So in this subsection, we discuss the influence of the squeezefilm damping on the dynamic behavior of variable geometry microbeam.We adopt the following function as squeeze-film damping force [25]: where  eff is the effective viscosity coefficient.
Here we also choose small voltage  dc = 10 V and  ac = 0.5 V.The response curves show hardening behavior.
Figure 15 shows the difference of dynamic behavior of the variable geometry microbeam between no gas damping and considering the squeeze-film damping.It can be found that the effect of the squeeze-film damping on the dynamic behavior of the microbeam is obvious.If the device is not working under vacuum, modeling should consider the influence of the damping.Meanwhile, by adjusting the initial gap between the microbeam and the substrate, we find that the influence of the initial gap on the squeeze-film damping is also obvious.When the initial gap increases, the electrostatic force decreases.If not, consider the effect of squeeze-film damping, and the amplitude of the dynamic response of the microbeam should decrease.However, as shown in Figure 15, under the circumstance of considering the squeeze-film damping, the amplitude of the dynamic response of the microbeam increases slightly.This is because augmenting the initial gap makes the decrease velocity of the squeeze-film damping quicker than that of the electrostatic force.So the trend shows like that in Figure 15.

Fringing Effects.
For a prismatic fixed-fixed microbeam, the fringing effect has always been neglected.But for the variable geometry microbeam, we should take the fringing effect into consideration because as the midsection of the microbeam becomes narrower, the influence of the fringing effect increases accordingly.This can be seen from Figure 16.The electric field energy leakage in the midsection of the variable geometry microbeam is larger than the counterpart prismatic beam and the fringing effect is more obvious.Thus, it is necessary to consider the influence of the fringing effect when calculating the dynamic characteristics of the variable geometry microbeam.Many different formulas for computing fringing fields have been proposed in [26,27] and we compared the difference of dynamic behavior of the microbeam among four typical fringing effect models.The revised formulas of electrostatic force  of different models are depicted in appendix.
Palmer's model is a classical one, which considered firstly the flux passing between the backsides of the plates.He adopted the Schwarz-Christoffel transformation method.But when applied to a narrow microbeam, the result of this model is poor.Fokkema's model took the effect of both the finite width and finite thickness of the microbeam into consideration.It can provide more accurate estimation for narrow microbeams.However, it also has limitations.When the gap size is larger than the beam width, the estimate is worse.Batra's model recalculated and optimized the parameters of empirical formula for the capacitance.It has a broader scope than Fokkema's model, especially for those significantly narrow microbeams.On the basis of Yang [28] and Palmer's research, Leus et al. proposed a new modified formula (A.4).This modified estimate is much more accurate than those in [28] and the relative error is very small compared with finite elements.
Figure 17 shows the response curves of forward and backward frequency sweeps.Among all these cases, the amplitude of the parallel-plate approximation (without regard to the fringing effects) is the maximum.Others, from high to low, are in turns Leus's model, Palmer's model, Batra's model, and Fokkema's model.As mentioned above, Batra and Fokkema's models are more suitable for narrow microbeam.However, the microbeam used in this work cannot be completely regarded as narrow microbeam because only the midsection is narrow enough.Meanwhile, the response results of Palmer and Leus's models are similar.Considering the calculation time, Palmer's model is selected as the compute model in this work.

Conclusions
Generalized differential quadrature method is an efficient method for the analysis of electrostatically actuated microbeam based devices.In this paper, a set of optimized variable geometry microbeams, accounting for the system nonlinearities due to midplane stretching, electrostatic force, and fringing effects, are analyzed to study their static and dynamic characteristics.We find that the more gradual change in width, the larger the maximum amplitude at resonance.Axial stress has prominent effect on the resonant frequency tunability characteristics.
The influence of different geometry shapes on the dynamic response is investigated.We find that the amplitude of the frequency response and the value of parameter  are positive correlation, but the discrepancy between the response curves is very small and with the increasing of , the discrepancy tends to be smaller.Augmenting  increases both the amplitude and resonant frequency of the microbeam, while the amplitude and the values of  are negative correlation.When  is small, the change of  has the similar trend with , but when the value of  is larger than 1.4, the trend becomes the opposite.In summary, , , and  should be changed when requiring an obvious frequency shift, while  can be used in fine tuning.In particular, when  is chosen as the variable, the designer should consider the special phenomenon mentioned above.
In the end, we discussed the influence of the squeeze-film damping and different fringing effect models on the dynamic behavior of the variable geometry microbeam.We found that, under the circumstance of considering the squeezefilm damping, the amplitude of the dynamic response of the microbeam increases slightly as the increase of the initial gap.

Figure 2 :
Figure 2: Static normalized midpoint deflection  max with DC voltage   .

Figure 4 :
Figure 4: Comparison of the frequency response curve: (a) in the hardening domain ( dc = 2 V,  ac = 0.1 V) and (b) in the softening domain ( dc = 3.5 V,  ac = 0.1 V).

Figure 11 :Figure 12 :Figure 13 :
Figure 11: Comparison of dynamic response curves of different midsection width .

Figure 14 :
Figure 14: Frequency response curve of different values of residual stress. = 0 means no residual stress is applied to the microbeam.

Figure 15 :
Figure 15: Frequency response curve of different initial gaps under squeeze-film damping.

Figure 16 :Figure 17 :
Figure 16: Electric field distribution and the corresponding contour of electric potential: (a) electric field distribution of prismatic beam; (b) electric field distribution of variable geometry microbeam; (c) contour of electric potential of prismatic beam; (d) contour of electric potential of variable geometry microbeam.

Table 1 :
Comparison of pull-in parameters of three available data sets.

Table 4 :
Dimension and physical parameters of a typical fixed-fixed microbeam.