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

Abstract

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.

1. 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–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.

2. 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 initial-boundary-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 𝑚th order partial derivative of 𝑓(𝑧,𝑡) with respect to 𝑧. For a function 𝑓(𝑧,𝑡), the differential quadrature approximation for the 𝑚th order derivative at the 𝑖th sampling point is given by 𝜕^𝑚𝜕𝑧^𝑚⎧⎪⎨⎪⎩𝑓(𝑧¹𝑓(𝑧₂,𝑡)𝑓(𝑧⋮_𝑁,𝑡),𝑡)⎫⎪⎬⎪⎭≅𝐷^{(𝑚)𝑖𝑗}⎧⎪⎨⎪⎩𝑓(𝑧¹𝑓(𝑧₂,𝑡)𝑓(𝑧⋮_𝑁,𝑡)for𝑖,𝑗=1,2,…𝑁,,𝑡)⎫⎪⎬⎪⎭(1) in which 𝑓(𝑧_𝑖,𝑡) is the functional value at the sample point 𝑧_𝑖, and 𝐷^{(𝑚)𝑖𝑗} are the differential quadrature weighting coefficients of the 𝑚th 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 𝐷^{(𝑚)𝑖𝑗}, that is, 𝑓(𝑧,𝑡)=^𝑁_𝑖=1𝑀(𝑧)(𝑧−𝑧_𝑖)𝑀₁(𝑧_𝑖)𝑓(𝑧^𝑖,𝑡),(2) where 𝑀(𝑧)=^𝑁_𝑗=1𝑀(𝑧−𝑧_𝑗),₁(𝑧_𝑖)=^𝑁_𝑗=_{1,𝑗≠𝑖}(𝑧_𝑖−𝑧_𝑗)for𝑖=1,2,…,𝑁.(3) Equation (2) is substituted into (1). The differential quadrature weighting coefficients are then given as 𝐷^(1)𝑖𝑗=𝑀¹(𝑧_𝑖)𝐷(𝑧_𝑖−𝑧_𝑗)𝑀₁(𝑧_𝑗)for𝑖,𝑗=1,2,…,𝑁,𝑖≠𝑗,_(1)𝑖𝑖=−^𝑁_𝑗=_{1,𝑗≠𝑖}𝐷^(1)𝑖𝑗for𝑖=1,2,…,𝑁.(4) Once the sampling points, such as 𝑧_𝑖 for 𝑖=1,2,…,𝑁, 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 𝐷𝐷^(2)𝑖𝑗=∑^𝑁_𝑘=1𝐷^(1)𝑖𝑘𝐷^(1)𝑘𝑗for𝑖,𝑗=1,2,…,𝑁,_(3)𝑖𝑗𝐷=∑^𝑁_𝑘=1𝐷^(1)𝑖𝑘𝐷^(2)𝑘𝑗for𝑖,𝑗=1,2,…,𝑁,_(4)𝑖𝑗𝐷=∑^𝑁_𝑘=1𝐷^(1)𝑖𝑘𝐷^(3)𝑘𝑗⋮for𝑖,𝑗=1,2,…,𝑁,_{(𝑚)𝑖𝑗}=∑^𝑁_𝑘=1𝐷^(1)𝑖𝑘𝐷^{(𝑚−1)𝑘𝑗}for𝑖,𝑗=1,2,…,𝑁.(5) 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 𝑧_𝑖=𝐿21−cos(𝑖𝑁−−1)1𝜋for𝑖=1,2,…,𝑁.(6) 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.

3. 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. 𝑡₀ specifies the thickness at the root of the actuator. 𝐿 is the length of the microactuator. 𝑞(𝑧) is the load. Load 𝑞(𝑧) acts on 𝑧=0∼𝐿 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 pull-in 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 𝑇=12^𝐿0𝜌𝐴𝜕𝑢𝜕𝑡²𝑑𝑧+12^𝐿0𝜌𝐴𝜕𝑣𝜕𝑡²+1𝑑𝑧2^𝐿0𝐺𝐽_𝑧𝜕𝜙𝜕𝑡²𝑑𝑧,(7) where 𝑢 is the displacement in the direction of the 𝑥-axis, 𝑣 is the displacement in the direction of the 𝑦-axis, 𝜙 is the twist angle in the direction of the 𝑧-axis, 𝐴 is the area of the cross-section of the microbeam, 𝐽_𝑧 is the polar moment in the direction of the 𝑧-axis, and 𝜌 is the density of the material of the actuator.

Figure 1 

Schematic view of an electrostatic fixed-fixed actuator.

While the external voltage 𝑒 is applied between the deformable beam and the fixed electrode, a position-dependent 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 𝑈=12^𝐿0𝐸𝐼_𝑦𝑦𝜕²𝑢𝜕𝑧²²+2𝐼_𝑥𝑦𝜕𝑢𝜕𝑧𝜕𝜕𝑣𝑧+𝐼^𝑥𝑥𝜕²𝑣+1𝜕𝑧²²𝑑𝑧2^𝐿0𝐺𝐽_𝑧𝜕𝜙𝜕𝑧²−1𝑑𝑧2^𝐿0𝑃𝜕𝑢𝜕𝑧²𝑑𝑧−12^𝐿0𝑃𝜕𝑣𝜕𝑧²𝑑𝑧,(8) where 𝐸 is Young’s modulus of the actuator, 𝐺 is the shear modulus, and 𝐼_𝑥𝑥, 𝐼_𝑦𝑦, and 𝐼_𝑥𝑦 are the moments of area. The load 𝑃 is the residual axial loading acting on the fixed end of the actuator. The value of 𝑃 is 𝜎𝑏₀𝑡₀. 𝜎 is the residual stress, and 𝑏₀ 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 𝑐_𝑢(𝜕𝑢/𝜕𝑡), 𝑐_𝑣(𝜕𝑣/𝜕𝑡), and 𝑐_𝜙(𝜕𝜙/𝜕𝑡) are assumed for resistance to the transverse velocity of the actuator. The damping forces 𝑐_𝑢𝑖(𝜕²/𝜕𝑧²)(𝐸𝐼(𝜕³𝑢/𝜕𝑡𝜕𝑧²)), 𝑐_𝑣𝑖(𝜕²/𝜕𝑧²)(𝐸𝐼(𝜕³𝑣/𝜕𝑡𝜕𝑧²)), and 𝑐_𝜙𝑖(𝜕/𝜕𝑧)(𝐺𝐽_𝑧(𝜕²𝜙/𝜕𝑡𝜕𝑧)) 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 𝛿𝑊 done by the bent actuator is 𝛿𝑊=−^𝐿0𝑐_𝑢𝜕𝑢𝜕𝑡𝛿𝑢𝑑𝑧−^𝐿0𝑐_𝑢𝑖𝜕²𝜕𝑧²𝐸𝐼𝜕³𝑢−𝜕𝑡𝜕𝑧²𝛿𝑢𝑑𝑧_𝐿0𝑐_𝑣𝜕𝑣𝜕𝑡𝛿𝑣𝑑𝑧−^𝐿0𝑐_𝑣𝑖𝜕²𝜕𝑧²𝐸𝐼𝜕³𝑣−𝜕𝑡𝜕𝑧²𝛿𝑣𝑑𝑧_𝐿0𝑐_𝜙𝜕𝜙𝜕𝑡𝛿𝜙𝑑𝑧+^𝐿0𝑐_𝜙𝑖𝜕𝜕𝑧𝐺𝐽^𝑧𝜕²𝜙+𝜕𝑡𝜕𝑧𝛿𝜙𝑑𝑧_𝐿0𝜀₀𝑏₀𝑒²2(𝑑+𝑆−𝛼sin(𝜋𝑧/𝐿)𝑡₀/2−𝑣)²−𝛿𝑣𝑑𝑧_𝐿0𝑞(𝑧)𝛿𝑣𝑑𝑧,(9) where 𝑒 is the applied voltage, 𝜀₀ is the dielectric constant of air, such as 𝜀₀=8.85×10⁻¹², 𝑏₀ is the width of the actuator, and 𝑑 is the initial gap as shown in Figure 1. The cross-section area of the actuator is 𝐴(𝑧)=𝑏₀𝑡₀(1+𝛼sin(𝜋𝑧/𝐿)), 𝛼 is the constant. 𝐼(𝑧) is the moment of inertia of the cross-sectional area of the actuator, which is 𝐼(𝑧)=𝐼₀(1+𝛼sin(𝜋𝑧/𝐿))³ and 𝐼₀=𝑏₀𝑡³⁰/12. The shape function 𝑆(𝑧) describes the shape of the curved electrode, and it is presented as 𝑆(𝑧)=𝛿_𝑒+𝛽sin(𝜋𝑧/𝐿). 𝛿_𝑒 is the fixed end gap distance of the curved electrode at 𝑧=0 and 𝑧=𝐿. The electrode shape is varied with the values of 𝛽 and 𝛿_𝑒. 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: ^𝑡2_𝑡1(𝛿𝑇−𝛿𝑈+𝛿𝑊)𝑑𝑡=0,(10) the dynamic deflection of a fixed-fixed micro-actuator can be expressed as the following nonlinear differential equation: 𝐸𝜕²𝐼_𝑦𝑦𝜕𝑧²𝜕²𝑢𝜕𝑧²+2𝐸𝜕𝐼^𝑦𝑦𝜕𝑧𝜕³𝑢𝜕𝑧³+𝐸𝐼_𝑦𝑦𝜕⁴𝑢𝜕𝑧⁴+𝐸𝜕²𝐼_𝑥𝑦𝜕𝑧²𝜕²𝑣𝜕𝑧²+2𝐸𝜕𝐼^𝑥𝑦𝜕𝑧𝜕³𝑣𝜕𝑧³+𝐸𝐼_𝑥𝑦𝜕⁴𝑣𝜕𝑧⁴+𝑃𝜕²𝑣𝜕𝑧²+𝑐_𝑢𝜕𝑢𝜕𝑡+𝑐^𝑢𝑖𝜕²𝜕𝑧²𝐸𝐼^𝑦𝑦𝜕³𝑢𝜕𝑧²𝜕𝑡+𝜌𝐴𝜕²𝑢𝜕𝑡²𝐸𝜕=0,₂𝐼_𝑥𝑥𝜕𝑧²𝜕²𝑣𝜕𝑧²+2𝐸𝜕𝐼^𝑥𝑥𝜕𝑧𝜕³𝑣𝜕𝑧³+𝐸𝐼_𝑥𝑥𝜕⁴𝑣𝜕𝑧⁴+𝐸𝜕²𝐼_𝑥𝑦𝜕𝑧²𝜕²𝑢𝜕𝑧²+2𝐸𝜕𝐼^𝑥𝑦𝜕𝑧𝜕³𝑢𝜕𝑧³+𝐸𝐼_𝑥𝑦𝜕⁴𝑢𝜕𝑧⁴+𝑃𝜕²𝑣𝜕𝑧²+𝑐_𝑣𝜕𝑣𝜕𝑡+𝑐^𝑣𝑖𝜕²𝜕𝑧²𝐸𝐼^𝑥𝑥𝜕³𝑣𝜕𝑧²𝜕𝑡+𝜌𝐴𝜕²𝑣𝜕𝑡²=𝜀⁰𝑏₀𝑒²2(𝑑+𝑆(𝑧)−𝛼sin(𝜋𝑧/𝐿)𝑡₀/2−𝑣(𝑧))²−𝜕−𝑞(𝑧)𝜕𝑧𝐺𝐽^𝑧𝜕𝜙𝜕𝑧+𝑐^𝜙𝜕𝜙𝜕𝑡−𝑐^𝜙𝑖𝜕𝜕𝑧𝐺𝐽^𝑧𝜕²𝜙+𝜌𝐽_𝑧𝜕𝑡𝜕𝑧𝜕²𝜙𝜕𝑡²=0,(11) where 𝜀₀ is the dielectric constant of air. The corresponding boundary conditions of the clamped-clamped micro-ctuator are 𝜕𝑢𝑢((00,,𝑡𝑡))=0,𝜕𝑢𝑢((𝜕𝐿𝐿𝑧,,𝑡𝑡))==00,,𝜕𝑣𝑣(𝜕(00𝑧,,𝑡𝑡))==00,,𝜕𝑣𝑣((𝜕𝐿𝐿𝑧,,𝑡𝑡))==00,,𝜙𝜙((𝜕𝐿0𝑧,,𝑡𝑡))===000,.,(12) 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 𝜕²𝑤𝜕𝑡²+[𝐶]𝜕𝜕𝑤𝑡+[𝐾]{𝑤}={𝐹}.(13) The displacement vector at the sample points is {𝑤}=⎧⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎩𝑢(𝑧¹𝑢(𝑧₂)𝑢(⋮𝑧_𝑁)𝑣(𝑧₁)𝑣(𝑧₂)𝑣(⋮𝑧_𝑁)𝜙(𝑧₁)𝜙(𝑧₂)𝜙(⋮𝑧_𝑁))⎫⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎭.(14)

The elements in the mass matrix are 𝑀𝑀_𝑖𝑖=0for𝑖=1,2,_𝑖𝑖𝑀=𝜌𝐴for𝑖=3,4,…,𝑁−2,_𝑖𝑖𝑀=0for𝑖=𝑁−1,𝑁,_𝑖𝑗𝑀=0for𝑖≠𝑗,𝑖=1,2,…,𝑁,𝑗=1,2,…,𝑁,_𝑖𝑗𝑀=0for𝑖=1,2,…,𝑁,𝑗=𝑁+1,𝑁+2,…,2𝑁,_𝑖𝑗𝑀=0for𝑖=1,2,…,𝑁,𝑗=2𝑁+1,2𝑁+2,…,3𝑁,_𝑖𝑗𝑀=0for𝑖=𝑁+1,𝑁+2,…,2𝑁,𝑗=1,2,…,𝑁,_𝑖𝑖𝑀=0for𝑖=𝑁+1,𝑁+2,_𝑖𝑖𝑀=𝜌𝐴for𝑖=𝑁+3,𝑁+4,…,2𝑁−2,_𝑖𝑖𝑀=0for𝑖=2𝑁−1,2𝑁,_𝑖𝑗𝑀=0for𝑖≠𝑗,𝑗𝑖==𝑁𝑁++11,,𝑁𝑁++22,,……,,22𝑁𝑁,_𝑖𝑗,𝑀=0f𝑗or=𝑖2=𝑁𝑁++1,12,𝑁𝑁++22,,……,,32𝑁𝑁,_𝑖𝑗,𝑀=0for𝑖=2𝑁+1,𝑗2𝑁=1+,22,,……,,𝑁3𝑁,_𝑖𝑗,𝑀=0for𝑖=𝑗=2𝑁𝑁++11,,2𝑁𝑁++22,,……,,23𝑁𝑁,_𝑖𝑖,𝑀=0for𝑖=2𝑁+1,_𝑖𝑖𝑀=𝜌𝐽_𝑧for𝑖=2𝑁+2,2𝑁+3,…,3𝑁−1,_𝑖𝑖𝑀=0for𝑖=3𝑁,_𝑖𝑗=0for𝑖≠𝑗,𝑗𝑖==22𝑁𝑁++11,,22𝑁𝑁++22,,……,,33𝑁𝑁.,(15) The elements in the damping matrix are 𝐶𝐶_𝑖𝑗=0for𝑖=1,2,𝑗=1,2,…,3𝑁,_𝑖𝑖=𝑐_𝑢+𝑐_𝑢𝑖𝜕²𝜕𝑧²(𝐸𝐼_𝑦𝑦)𝐷^(2)𝑖𝑖+2𝑐_𝑢𝑖𝜕𝜕𝑧(𝐸𝐼^𝑦𝑦)𝐷^(3)𝑖𝑖𝐶+𝑐_𝑢𝑖𝐸𝐼_𝑦𝑦𝐷^(4)𝑖𝑖for𝑖=3,4,…,𝑁−2,_𝑖𝑗=𝑐_𝑢𝑖𝜕²𝜕𝑧²(𝐸𝐼_𝑦𝑦)𝐷^(2)𝑖𝑗+2𝑐_𝑢𝑖𝜕𝜕𝑧(𝐸𝐼^𝑦𝑦)𝐷^(3)𝑖𝑗+𝑐_𝑢𝑖𝐸𝐼_𝑦𝑦𝐷^(4)𝑖𝑗𝐶for𝑖≠𝑗,𝑖=3,4,…,𝑁−2,𝑗=1,2,…,𝑁,_𝑖𝑗𝐶=0for𝑗𝑖==3𝑁,4+,…1,𝑁,𝑁+−23,…,𝑁,−3𝑁2,_𝑖𝑗,𝐶=0for𝑖=𝑁−1,𝑁,𝑗=1,2,…,3𝑁,_𝑖𝑗𝐶=0for𝑖=𝑁+1,𝑁+2,…,2𝑁,𝑗=1,2,…,𝑁,_𝑖𝑗𝐶=0for𝑖=𝑁+1,𝑁+2,𝑗=𝑁+1,𝑁+2,…,2𝑁,_𝑖𝑖=𝑐_𝑣+𝑐_𝑣𝑖𝜕²𝜕𝑧²(𝐸𝐼_𝑥𝑥)𝐷^(2)𝑖𝑖+2𝑐_𝑣𝑖𝜕𝜕𝑧(𝐸𝐼^𝑥𝑥)𝐷^(3)𝑖𝑖+𝑐_𝑣𝑖𝐸𝐼_𝑥𝑥𝐷^(4)𝑖𝑖𝐶for𝑖=𝑁+3,𝑁+4,…,2𝑁−2,_𝑖𝑗=𝑐_𝑣𝑖𝜕²𝜕𝑧²(𝐸𝐼_𝑥𝑥)𝐷^(2)𝑖𝑗+2𝑐_𝑣𝑖𝜕𝜕𝑧(𝐸𝐼^𝑥𝑥)𝐷^(3)𝑖𝑗+𝑐_𝑣𝑖𝐸𝐼_𝑥𝑥𝐷^(4)𝑖𝑗𝐶for𝑖≠𝑗,𝑖=𝑁+3,𝑁+𝑗4,=…𝑁,2+𝑁1−,𝑁2+,2,…,2𝑁,_𝑖𝑗𝐶=0for𝑖=2𝑁−1,2𝑁,𝑗=𝑁+1,𝑁+2,…,2𝑁,_𝑖𝑗𝐶=0f𝑗or=𝑖2=𝑁𝑁++1,12,𝑁𝑁++22,,……,,32𝑁𝑁,_𝑖𝑗,𝐶=0for𝑖=2𝑁+1,𝑗=1,2,…,3𝑁,_𝑖𝑗𝐶=0for𝑖=2𝑁+2,2𝑁𝑗+=31,,…2,…,3𝑁,2−𝑁,_𝑖𝑖1,=𝑐_𝜙−𝑐_𝜙𝑖𝜕𝜕𝑧(𝐺𝐽^𝑧)𝐷^(1)𝑖𝑖−𝑐_𝜙𝑖(𝐺𝐽_𝑧)𝐷^(2)𝑖𝑖𝐶for𝑖=2𝑁+2,2𝑁+3,…,3𝑁−1,_𝑖𝑗=−𝑐_𝜙𝑖𝜕𝜕𝑧(𝐺𝐽^𝑧)𝐷^(1)𝑖𝑗−𝑐_𝜙𝑖(𝐺𝐽_𝑧)𝐷^(2)𝑖𝑗𝐶for𝑖≠𝑗,𝑖=𝑗2=𝑁2+𝑁2+,21𝑁,2+𝑁3+,…2,,…3𝑁,3−𝑁1,_𝑖𝑗,=0for𝑖=3𝑁,𝑗=1,2,…,3𝑁.(16) The elements in the stiffness matrices are 𝐾𝐾₁₁=1,_1𝑗𝐾=0for𝑗=2,3,…,3𝑁,_2𝑗=𝐷^(1)1𝑗𝐾𝐿for𝑗=1,2,…,𝑁,_2𝑗𝐾=0for𝑗=𝑁+1,𝑁+2,…,3𝑁,_𝑖𝑗=𝐸𝜕²𝐼_𝑦𝑦𝜕𝑧²|||_{𝑧=𝑧𝑖}𝐷^(2)𝑖𝑗+2𝐸𝜕𝐼^𝑦𝑦𝜕𝑧|||_{𝑧=𝑧𝑖}𝐷^(3)𝑖𝑗𝐾+𝐸𝐼_𝑦𝑦𝐷^(4)𝑖𝑗+𝑃𝐷^(2)𝑖𝑗for𝑖=𝑗3=,41,,…2,,…𝑁,−𝑁2,_𝑖𝑗,=𝐸𝜕²𝐼_𝑥𝑦𝜕𝑧²|||_{𝑧=𝑧𝑖}𝐷^{(2)𝑖,𝑗−𝑁}+2𝐸𝜕𝐼^𝑥𝑦𝜕𝑧|||_{𝑧=𝑧𝑖}𝐷^{(3)𝑖,𝑗−𝑁}𝐾+𝐸𝐼_𝑥𝑦𝐷𝑗^{(4)𝑖,𝑗−𝑁}=𝑁fo+r𝑖1,=𝑁3+,42,…,…,𝑁,2𝑁−,_{𝑁−1,𝑗}2,𝐾=𝐷^{(1)𝑁,𝑗−𝑁}for𝑗=1,2,…,𝑁,_𝑁𝑗=0fo𝐾r𝑗=1,2,…,𝑁−1,_𝑁𝑁𝐾=1,_𝑁𝑗𝐾=0for𝑗=𝑁+1,𝑁+2,…,3𝑁,_𝑖𝑗𝐾𝑗==02𝑁fo+r𝑖1=,2𝑁3,4+,…2,…,𝑁,3−𝑁2,_𝑖𝑗,=0for𝐾𝑖=𝑁+1,𝑗=1,2,…,𝑁,_{𝑁+1,𝑁+1}𝐾=1,_𝑖𝑗=𝐾0for𝑖=𝑁+1,𝑖=𝑁+2,𝑁+3,…,3𝑁,_𝑖𝑗𝐾=0for𝑖=𝑁+2,𝑗=1,2,…,𝑁,_𝑖𝑗𝐾=𝐷^{(1)1,𝑗−𝑁}for𝑖=𝑁+2,𝑗=𝑁+1,𝑁+2,…,2𝑁,_𝑖𝑗𝐾=0for𝑖=𝑁+2,𝑖=2𝑁+1,2𝑁+2,…,3𝑁,_𝑖𝑗=𝐸𝜕²𝐼_𝑥𝑦𝜕𝑧²|||_{𝑧=𝑧𝑖}𝐷^{(2)𝑖−𝑁,𝑗}+2𝐸𝜕𝐼^𝑥𝑦𝜕𝑧|||_{𝑧=𝑧𝑖}𝐷^{(3)𝑖−𝑁,𝑗}+𝐸𝐼_𝑦𝑦𝐷^(4)𝑖𝑖𝐾for𝑖=𝑁+3,𝑁+4,…,2𝑁−2,𝑗=1,2,…,𝑁,_𝑖𝑗=𝐸𝜕²𝐼_𝑥𝑥𝜕𝑧²|||_{𝑧=𝑧𝑖}𝐷^{(2)𝑖−𝑁,𝑗−𝑁}+2𝐸𝜕𝐼^𝑥𝑥𝜕𝑧|||_{𝑧=𝑧𝑖}𝐷^{(3)𝑖−𝑁,𝑗−𝑁}+𝐼_𝑥𝑥𝐷^{(4)𝑖−𝑁,𝑗−𝑁}𝐾for𝑖=𝑁+3,𝑖𝑁=+𝑁4,+…1,,𝑁2𝑁+−2,2…,,2𝑁,_𝑖𝑗𝐾=0for𝑖=2𝑁−1,𝑗=1,2,…,𝑁,_𝑖𝑗𝐾=𝐷^{(1)1,𝑗−𝑁}for𝑖=2𝑁−1,𝑗=𝑁+1,𝑁+2,…,2𝑁,_𝑖𝑗=𝐾0for𝑖=2𝑁−1,𝑗=2𝑁+1,2𝑁+2,…,3𝑁,_𝑖𝑗𝐾=0for𝑖=2𝑁+1,𝑗=1,2,…,2𝑁,_𝑖𝑗𝐾=1for𝑖=2𝑁+1,𝑗=2𝑁+1,_𝑖𝑗𝐾=0for𝑖=2𝑁+1,𝑗=2𝑁+2,…,3𝑁,_𝑖𝑗𝐾=0for𝑖=2𝑁+2𝑗,2=𝑁1+,23,…,…,2,3𝑁𝑁,_𝑖𝑗,=−𝐺𝜕𝐽^𝑧𝜕𝑧|||_{𝑧=𝑧𝑖}𝐷^{(1)𝑖−2𝑁,𝑗−2𝑁}−𝐺𝐽_𝑧𝐷^{(2)𝑖−2𝑁,𝑗−2𝑁}𝐾for𝑖=𝑗2𝑁=2+𝑁2,+2𝑁1,2+𝑁3,+…2,,…3𝑁,3−𝑁1,_𝑖𝑗,=𝐾0for𝑖=3𝑁,𝑗=1,2,…,3𝑁−1,_𝑖𝑗𝐹=1for𝑖=3𝑁,𝑗=3𝑁,_𝑖𝐹=0for𝑖=1,2,…,𝑁+2,_𝑖=𝜀⁰𝑏𝑒²2(𝑑+𝑆(𝑧)−𝛼sin(𝜋𝑧/𝐿)𝑡₀/2−𝑣(𝑧))²𝐹for𝑖=𝑁+3,𝑁+4,…,2𝑁−−𝑞(2𝑧,_𝑖)=0for𝑖=2𝑁−1,2𝑁,…,3𝑁.(17) 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-by-step 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.

4. Numerical Results and Discussion

The microactuator is fabricated from polysilicon material. The geometric parameters and the material of the microactuator are 𝐸=150GPa, 𝛿_max=30𝜇m, 𝛼=0, 𝛽=0, 𝑐_𝑢𝑖=0, 𝑐_𝑣𝑖=0, 𝑐_𝜙𝑖=0, 𝑐_𝑢=0, 𝑐_𝑣=0, 𝑐_𝜙=0, 𝑏₀=5𝜇m, 𝑡₀=2𝜇m, 𝑑=2𝜇m, and 𝐿=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 fixed-fixed actuator for various residual stresses. The value of applied voltageis 620V. The nonlinear dynamic equation formed by the differential quadrature method is solved by the Wilson-𝜃 integration method, with 𝜃=1.4 and Δ𝑡=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 fixed-fixed 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.

Figure 2 

Deflections of an electrostatic fixed-fixed actuator for various applied voltages.

Figure 3 

Frequencies of an electrostatic fixed-fixed actuator for various beam lengths.

Figure 4 

Deflections near the middle of an electrostatic fixed-fixed actuator for various residual stresses.

Figure 5 

Stresses near the middle of an electrostatic fixed-fixed actuator for various residual stresses.

Figure 6 

Stresses near the root of an electrostatic fixed-fixed actuator for various residual stresses.

5. 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 pull-in phenomena have also been investigated by employing the proposed differential quadrature method algorithm.