The Effect of Residual Stress on the Electromechanical Behavior of Electrostatic Microactuators

This work simulates the nonlinear electromechanical behavior of different electrostatic microactuators. It applies the differential quadrature method, Hamilton’s principle, and Wilson-θ integration method to derive the equations of motion of electrostatic microactuators and find a solution to these equations. Nonlinear equation difficulties are overcome by using the differential quadrature method. The stresses of electrostatic actuators are determined, and the residual stress effects of electrostatic microactuators are simulated.


INTRODUCTION
Osterberg et al. [1] analyzed electrostatically deformed diaphragms using a one-dimensional numerical model and a three-dimensional model.Osterberg and Senturia [2] showed that the sharp instability phenomena of electrostatic pull-in behavior of cantilever beam and fixed-fixed beam actuators can be adopted to extract the material properties of microelectromechanical system.Elwenspoek et al. [3] studied the dynamic behavior of active joints for various electrostatic actuator designs.Hirai et al. [4][5][6] presented the deflection characteristics of electrostatic actuators with modified electrode and cantilever shapes.Wang [7] applied a feedback control for suppressing the vibration of actuator beams in an electrostatic actuator.Shi et al. [8] combined an exterior boundary element method for electrostatics and a finite element method for elasticity to evaluate the coupling effect between the electrostatic force and the elastic deformation.Gretillat et al. [9] employed three-dimensional finite element programs to simulate the dynamics of a nonlinear actuator, considering the effect of squeeze-film damping.Hung and Senturia [10] proposed leveraged bending and strain-stiffening methods to enlarge the limit of travel distance before pull-in of electrostatic actuators.This work will analyze the nonlinear pull-in behaviors of different types of microactuator with various residual stresses using the differential quadrature method.The Chebyshev-Gauss-Lobatto point distribution on each actuator will be used.The integrity and computational accuracy of the differential quadrature method in solving this problem will be evaluated through a range of case studies.The dynamic equations of the cantilever microactuator are derived using the differential quadrature method.The equations describing the residual vibrations of the microelectrostatic actuators are derived in this paper.The differential quadrature method is used to produce the electrostatic field equations in matrix form.

THE DIFFERENTIAL QUADRATURE METHOD
This paper employs the differential quadrature method, with its easy-to-use and meshless technique, to analyze the nonlinear deflection behaviors of different types of microactuator with different residual stresses.There are a number of solution techniques for complicated beam problems, such as the Rayleigh-Ritz method, the analytical method, the Galerkin method, the finite element method, and the boundary element method.The differential quadrature method has been extensively used to solve a variety of problems in different fields of science and engineering with no need of energy formulation.The differential quadrature method has been shown to be a powerful contender in solving initial and boundary value problems and has, thus, become an alternative to the previous methods.Jang et al. [11] proposed the δ method, in which the boundary points are selected at a small distance from each other.The δ technique can be applied to the double boundary conditions of plate and beam problems.The accuracy of the solution depends on a sufficiently small δ.The boundary points are chosen at a small distance δ.The δ technique can be applied to the double boundary conditions of plate and beam problems.The use of δ at the boundary makes the matrix ill conditioned [11].Wang and Bert [12] considered the boundary conditions in finding the differential quadrature weighting coefficients.Malik and Bert [13] solved the problem of the free vibration of the plates and showed that the boundary conditions can be built into the differential quadrature weighting coefficients.In their formulation, the multiple boundary conditions are directly applied to the differential quadrature weighting coefficients, and thus it is not necessary to select a nearby point.In other words, the accuracy of the calculated results will be independent of the value of the δ-interval.The differential quadrature weighting coefficients can be obtained by multiplying the inverse matrix [13].Sherbourne and Pandey [14] solved buckling problems using the differential quadrature method.From the foregoing discussion, over the last two decades, the differential quadrature method has been applied extensively as an effective means of solving a range of problems in various fields of science and engineering.Quan and Chang [15,16] derived the weighting coefficients in a more explicit way.Feng and Bert [17] analyzed the flexural vibration analysis of a geometrically nonlinear beam using the quadrature method.Chen and Zhong [18] presented the study on the nonlinear computations of the differential quadrature method and differential cubature method.Tomasiello [19] applied the differential quadrature method to initialboundary-value problems.Wang et al. [20] presented the free vibration analysis of circular annular plates with nonuniform thickness by the differential quadrature method.Wang and Gu [21] presented the static analysis of frame structures by the differential quadrature element method.Liew et al. [22,23] presented the differential quadrature method for Mindlin plates on Winkler foundations and thick symmetric cross-ply laminates with first-order shear flexibility.Du et al. [24] presented the application of a generalized differential quadrature method to structural problems.Mirfakhraei and Redekop [25] solved the buckling of circular cylindrical shells using the differential quadrature method.Moradi and Taheri [26] presented the delamination buckling analysis of general laminated composite beams using the differential quadrature method.De Rosa and Franciosi [27] introduced the exact and approximate dynamic analysis of circular arches using the differential quadrature method.Sun and Zhu [28] used the upwind local differential quadrature method for solving incompressible viscous flow.Gu and Wang [29] presented the free vibration analysis of circular plates with stepped thickness over a concentric region by the differential quadrature method.Du et al. [30] presented the generalized differential quadrature method for buckling analysis.Han and Liew [31] analyzed axisymmetric free vibration of thick annular plates.Tanaka and Chen [32] applied a dual reciprocity boundaryelement method to transient elastodynamic problems using the differential quadrature method.Chen et al. [33] solved the high-accuracy plane stress and plate elements by the quadrature element method.The essence of the differential quadrature method is that the derivative of a function at a sample point can be approximated as a weighted linear summation of the functional values at all of the sampling points in the domain.Using this approximation, the differential equation is then reduced to a set of algebraic equations.The effects of position-dependent electrostatic force and axial residual stress have all been considered in the proposed models.While the efficiency and accuracy of the Rayleigh-Ritz method depend on the number and accuracy of the selected comparison functions; the differential quadrature method does not have this difficulty of selecting the appropriate comparison functions.The differential quadrature method approximates the mth order partial derivative of f (z, t) with respect to z.For a function f (z, t), the differential quadrature approximation for the mth order derivative at the ith sampling point is given by in which f (z i , t) is the functional value at the sample point z i , and D (m) i j are the differential quadrature weighting coefficients of the mth order differentiation attached to these functional values.Quan and Chang [15,16] introduced a Lagrangian interpolation polynomial to overcome the numerical ill conditioning in determining the differential quadrature weighting coefficients D (m)  i j , that is, where ( Equation ( 2) is substituted into (1).The differential quadrature weighting coefficients are then given as for i, j = 1, 2, . . ., N, i / = j, D (1)   ii = − N j=1, j / = i D (1)   i j Once the sampling points, such as z i for i = 1, 2, . . ., N, are selected, the coefficients of the differential quadrature weighting matrix can be obtained from (4).Higher order derivatives of the differential quadrature weighting coefficients can also be directly calculated by matrix multiplication [34], which can be expressed as D (1)  ik D (1)   k j for i, j = 1, 2, . . ., N, D (3)   i j = N k=1 D (1)  ik D (2)   k j for i, j = 1, 2, . . ., N, D (4)   i j = N k=1 D (1)  ik D (3)   k j for i, j = 1, 2, . . ., N, . . .
k j for i, j = 1, 2, . . ., N. ( There are many computational methods available for dynamic analysis.In this paper, the residual vibrations of the microelectrostatic actuators are investigated using the differential quadrature method.The most convenient approach to solving a beam structure problem is to uniformly space out the sample points.The selection of sample points is important for the accuracy of the differential quadrature method solution, but inaccurate results have been obtained when using this uniform distribution.A nonuniform sample point distribution, such as the Chebyshev-Gauss-Lobatto distribution [34], improves the accuracy of the calculation.The integrity and computational efficiency of the differential quadrature method in solving this problem will be demonstrated using a set of case studies.However, an alternative efficient technique is still sought.In this study, the unequally spaced sample points of each beam using the Chebyshev-Gauss-Lobatto distribution are chosen as The differential quadrature method has been shown to be a powerful candidate for solving initial and boundary value problems, and has thus become an alternative to other methods.The efficiency and the accuracy of Rayleigh-Ritz method depend on the number and accuracy of the selected comparison functions, whereas the differential quadrature method does not have such a difficulty.Like that of any polynomial approach, the accuracy of the solution using this method is improved by increasing the number of sample points.The differential quadrature method uses high-order element level, where the finite element method approximates a function using low-order polynomials.

DYNAMIC BEHAVIOR OF MICROACTUATORS
A shaped microbeam with a curved electrode is shown in Figure 1.The figure depicts the geometry of a tapered electrostatic microactuator.t 0 specifies the thickness at the root of the actuator.L is the length of the microactuator.q(z) is the load.Load q(z) acts on z = 0∼L in the beam.As a driving voltage is applied between the fixed-fixed microbeam and the electrode, a position-dependent electrostatic pressure is distributed to deform the microbeam toward the curved electrode.The gap between the shaped microbeam and the curved electrode determines the distribution of the electrostatic pressure.To prevent a short circuit after pullin contact, an isolated layer or other structure is required.
The force pulls the microbeam toward the shaped electrode.Different electrode shapes have been proposed to improve the electrostatic force distribution and the deformed shape of the actuator.The kinetic energy of the microactuator is where u is the displacement in the direction of the x-axis, v is the displacement in the direction of the y-axis, φ is the twist angle in the direction of the z-axis, A is the area of the cross-section of the microbeam, J z is the polar moment in the direction of the z-axis, and ρ is the density of the material of the actuator.While the external voltage e is applied between the deformable beam and the fixed electrode, a positiondependent electrostatic pressure is created to pull the deformable beam toward the ground electrode.This electrostatic pressure is approximately proportional to the inverse of the square of the gap between them.When the voltage reaches the critical voltage, the fixed-fixed beam will be pulled toward the electrode suddenly.The electric fringing effects are ignored in the following analyses.The strain energy of the microactuator can be approximated as where E is Young's modulus of the actuator, G is the shear modulus, and I xx , I yy , and I xy are the moments of area.The load P is the residual axial loading acting on the fixed end of the actuator.The value of P is σb 0 t 0 .σ is the residual stress, and b 0 is the beam width.Because of the coupling between the mechanical and electrostatic effects, the behavior of the electrostatic actuator appears more complicated than elastic behavior.The external damping presents a viscous resistance to transverse displacement of the actuator, and the internal damping provides a viscous resistance to straining of the microactuator material.The damping forces c u (∂u/∂t), c v (∂v/∂t), and c φ (∂φ/∂t) are assumed for resistance to the transverse velocity of the actuator.The damping forces ), and c φi (∂/∂z)(GJ z (∂ 2 φ/∂t∂z)) are assumed for the resistance to the strain velocity of the microactuator.Considering the electrostatic force and the internal and external damping effects in the actuator, the virtual work δW done by the bent actuator is where e is the applied voltage, ε 0 is the dielectric constant of air, such as ε 0 = 8.85 × 10 −12 , b 0 is the width of the actuator, and d is the initial gap as shown in Figure 1.The cross-section area of the actuator is A(z) = b 0 t 0 (1 + α sin(πz/L)), α is the constant.I(z) is the moment of inertia of the cross-sectional area of the actuator, which is I(z) = I 0 (1 + α sin(πz/L)) 3 and I 0 = b 0 t 3 0 /12.The shape function S(z) describes the shape of the curved electrode, and it is presented as S(z) = δ e + β sin(πz/L).δ e is the fixed end gap distance of the curved electrode at z = 0 and z = L.The electrode shape is varied with the values of β and δ e .However, due to the difficulty of no linearity between the actuator deflection and the electrostatic force, this residual vibration phenomenon has been studied in only a very few papers, as has the effect of electrode shape on the residual response.Substituting (7), (8), and (9) into Hamilton's equation: the dynamic deflection of a fixed-fixed micro-actuator can be expressed as the following nonlinear differential equation: where ε 0 is the dielectric constant of air.The corresponding boundary conditions of the clamped-clamped micro-ctuator are Equation ( 1) is substituted into ( 11)-( 12) by employing the differential quadrature method.The equations of motion of the microactuator can be discretized in matrix form with respect to the sample points as The displacement vector at the sample points is φ(z 2 ) . . .
The elements in the mass matrix are The elements in the damping matrix are ∂z 2 (EI yy )D (2)  ii + 2c ui ∂ ∂z (EI yy )D (3)   ii + c ui EI yy D (4)   ii (2)  i j + 2c ui ∂ ∂z (EI yy )D (3)  i j + c ui EI yy D (4)   i j ii + c vi EI xx D (4)   ii (1)  ii − c φi (GJ z )D (2)   ii The elements in the stiffness matrices are ∂I yy ∂z z=zi D (3)   i j + EI yy D (4)  i j + PD (2)   i j ∂I xy ∂z z=zi D (3)   i, j−N + EI xy D (4)   i, j−N for i = 3, 4, . . ., N − 2, N, j−N for j = 1, 2, . . ., N, K N j = 0 for j = 1, 2, . . ., N − 1, ∂I xy ∂z z=zi D (3)  i−N, j + EI yy D (4)   ii ∂I xx ∂z z=zi D (3)   i−N, j−N + I xx D (4)   i−N, j−N i−2N, j−2N − GJ z D (2)   i−2N, j−2N for i = 2N + 2, 2N + 3, . . ., 3N − 1, The dynamic responses of the microactuator are solved using the Wilson-θ integration method in this paper.The Wilson-θ integration method is an effective implicit time integration procedure for dynamic problems.It is a step-bystep integration method that assumes that the acceleration terms vary linearly between consecutive sampling instants.An electrostatic force pulls the cantilever actuator toward the curved electrode.The electrostatic force is generated by the difference between voltage applied to the curved electrode and that applied to the actuator.This electrostatic pressure is approximately proportional to the inverse of the square of the gap between them.When the voltage exceeds the critical voltage, the fixed-fixed beam is suddenly pulled into the electrode.

NUMERICAL RESULTS AND DISCUSSION
The microactuator is fabricated from polysilicon material.The geometric parameters and the material of the microactuator are E = 150 GPa, δ max = 30 μm, α = 0, β = 0, c ui = 0, c vi = 0, c φi = 0, c u = 0, c v = 0, c φ = 0, b 0 = 5 μm, t 0 = 2 μm, d = 2 μm, and L = 500 μm. Figure 2 shows the deflections of the microbeam with different applied voltages.The results indicate that the results calculated from the proposed differential quadrature method agree very well with the results found using the finite element method.Figure 3 shows the frequencies of an electrostatic fixed-fixed actuator for various lengths of the beam.Again, the results found using the differential quadrature method are similar to the results found using the finite element method.Figure 4 plots the deflections near the middle of an electrostatic fixedfixed actuator for various residual stresses.The value of applied voltageis 620 V.The nonlinear dynamic equation formed by the differential quadrature method is solved by  the Wilson-θ integration method, with θ = 1.4 and Δt = 0.003 millisecond.A number of papers state that Wilson-θ integration method is unconditionally stable with a factor of θ ≥ 1.37 [35,36].The calculated results show that higher residual stresses produce smaller deflections near the middle of an electrostatic fixed-fixed actuator.Figure 5 shows the stresses near the middle of an electrostatic fixed-fixed actuator for various residual stresses.Numerical results in this example show that the residual stresses can significantly affect the dynamic behavior of the actuator system, showing that higher residual stresses produce larger stresses near the middle of an electrostatic fixed-fixed actuator.Figure 6 shows the stress near the root of an electrostatic fixed-fixed actuator for various residual stresses.Results indicate that residual stress is a very sensitive parameter for the residual vibration of the microactuator.Numerical results in this example show that the driving voltage can affect the electromechanical behavior of the actuator system significantly.Calculated results also display that the higher residual stresses introduce the larger stresses near the root of an electrostatic fixedfixed actuator.Residual axial loading should be considered in the design.Numerical results indicate that the differential quadrature method is a feasible and efficient method to analyze the nonlinear pull-in behavior of a fixed-fixed type of electrostatic microbeam.

CONCLUSIONS
The differential quadrature method is highly suited to designing or analyzing an electrostatic microactuator.The simplicity of this formulation makes it a strong candidate for modeling applications that are more complicated.The effects of residual stresses of microactuators on the nonlinear pullin phenomena have also been investigated by employing the proposed differential quadrature method algorithm.

Figure 1 :
Figure 1: Schematic view of an electrostatic fixed-fixed actuator.

Figure 5 :Figure 6 :
Figure 5: Stresses near the middle of an electrostatic fixed-fixed actuator for various residual stresses.