# Dynamic Analysis of Electrostatic Microactuators Using the Differential Quadrature Method

Authors: Ming-Hung Hsu
Journal: Advances in Acoustics and Vibration
 (2011)

---

## Abstract

This work studies the dynamic behavior of electrostatic actuators using finite-element package software (FEMLAB) and differential quadrature method. The differential quadrature technique is used to transform partial differential equations into a discrete eigenvalue problem. Numerical results indicate that length, width, and thickness significantly impact the frequencies of the electrostatic actuators. The thickness could not affect markedly the electrostatic actuator capacities. The effects of varying actuator length, width, and thickness on the dynamic behavior and actuator capacities in electrostatic actuator systems are investigated. The differential quadrature method is an efficient differential equation solver.


---

## Body

## 1. Introduction

Plate-type electrostatic actuators are widely applied in microelectromechanical systems. Microelectrostatic actuator devices have a high operating frequency, low-power consumption and can replace many passive components. Mehdaoui et al. [1] presented the vertical cointegration of AlSi MEMS tunable capacitors and Cu inductors for tunable LC blocks. Etxeberria and Gracia [2] presented tunable MEMS volume capacitors for high-voltage applications. Liu et al. [3] presented that actuation by electrostatic repulsion is produced by nonvolatile charge injection. Gallant and Wood [4] investigated how fabrication techniques affect the performance of widely tunable micromachined capacitors. Borwick III et al. [5] analyzed a high Q, large tuning range MEMS capacitor for RF filter systems. Harsh et al. [6] studied the realization and design considerations of a flip-chip integrated MEMS tunable capacitor. Petersen [7] first described the nonlinear pull-in behavior of an electrostatic microactuator. Osterberg et al. [8] proposed different numerical models for analyzing electrostatically deformed diaphragms. Their results revealed that the electrostatic deformation calculated using the one-dimensional model is close to that obtained using a three-dimensional model. Various models, including the lumped parallel-plate spring model, the one-dimensional numerical model, and the finite-element model that incorporates a three-dimensional simulation, were proposed to calculate the pull-in behaviors of various fixed-fixed DMD structures and pressure sensors [8]. Gilbert et al. [9] analyzed the three-dimensional coupled electromechanics of MEMS using a CoSolve-EM simulation algorithm. Elwenspoek et al. [10] studied the dynamic behavior of active joints for various electrostatic actuator designs. Shi et al. [11] presented a combination of an exterior boundary element method for electrostatics combined with a finite-element method for elasticity to evaluate the effect of coupling between the electrostatic force and the elastic deformation. Osterberg and Senturia [12] applied the sharp instability phenomena of electrostatic pull-in behaviors for cantilever beam and fixed-fixed beam actuators to elicit the material characteristics of MEMS. Gretillat et al. [13] employed three-dimensional MEMCAD and FEM programs to simulate the dynamics of a nonlinear actuator, considering the effect of squeeze-film damping. Hung and Senturia [14] developed leveraged bending and strain-stiffening methods to increase the limiting travel distance before pull-in of electrostatic actuators. Chan et al. [15] measured the pull-in voltage and capacitance-voltage and performed 2D simulations that included the electrical effects of fringing fields and finite-beam thickness to determine the material properties of electrostatic microactuators. Li and Aluru [16] developed a mixed-regime approach for combining linear and nonlinear theories to analyze large MEMS deformations at large applied voltages. Chyuan et al. [17–19] established the validity and accuracy of the dual boundary element method and applied it to study the effect of gap size variation for the levitation of MEMS comb drive. Lai and Chen [20] studied the influence of the holes in the membrane structures of radiofrequency MEMS switches. Qiao et al. [21] presented the suspension beam called two-beam to achieve parallel-plate actuator with extended working range, but without penalties of complex control circuit and large actuation voltage. In this work, the finite-element method and the differential quadrature method are employed to analyze how actuator length, width, and thickness affect dynamic behavior and capacitances in electrostatic actuator systems.


## 2. Capacitance of Electrostatic Micro Actuators

Figure1 depicts the geometry of an electrostatic actuator where h is the thickness of the movable plate [22]. The length and width of the movable plate are a and b, respectively. The microactuator design is based upon the deformation of a movable mechanical structure by electrostatic forces using a fixed electrode. A voltage applied across the gap creates the electrostatic force. Electrostatic actuator devices are fabricated from polysilicon by surface micromaching techniques. In this work, the finite-element technique is applied to find the capacitance of the electrostatic microactuator. The finite-element method is used to find approximate solutions for partial differential equations and integral equations. The convergence and accuracy of the finite-element solution depends on the differential equation, integral form, and element used. The energy stored in the microdevice is expressed as follows [22–24]:
(1)We=12∫vρvVdv,
where ρv is the volume charge density, v is the volume, and V is the applied voltage:
(2)ρv=∇⃗⋅De⃗,
where De⃗ is the vector of electric flux density,
(3)C=2WeV2,
where C is the capacitance. The capacitance models are extracted using three-dimensional simulations. Figure 2 shows the capacitances of the electrostatic microactuator with various lengths. The width, thickness, and gap distance of the electrostatic actuator are 176, 8, and 38 micrometers, respectively. The computational solution is acquired by the finite-element scheme. The FEMLAB finite-element package is used to analyze the model. In previous cases, 62463 solid elements are used to calculate the capacitance of the actuator and provide convergent results. The effect of the actuator length is more pronounced at higher capacity. Figure 3 shows the capacitances of the electrostatic microactuator for varying widths. The length, thickness, and gap distance of the electrostatic actuator are 262, 8, and 38 micrometers, respectively. Numerical results reveal that the capacity of the actuator increases as the width of the actuator increases. Figure 4 shows the capacitances of electrostatic microactuator for various gap distances. The width, length, and thickness of the electrostatic actuator are 176, 262, and 8 micrometers, respectively. Numerical results reveal that the capacity of the actuator increases as the gap distance of the actuator falls. Figure 5 shows the capacitances of the electrostatic microactuator at various thicknesses. The width, length, and gap distance of the electrostatic actuator are 176, 262, and 38 micrometers, respectively. Numerical results show that the thickness does not significantly affect the capacitances of the electrostatic microactuator. The evaluation illustrates capacitance problems.Figure 1
Schematic view of electrostatic actuator.Figure 2
Capacities of microelectrostatic actuator for various lengths.Figure 3
Capacities of microelectrostatic actuator for various widths.Figure 4
Capacities of microelectrostatic actuator for various gap distances.Figure 5
Capacities of microelectrostatic actuator for various thicknesses.


## 3. Vibration Analysis of Electrostatic MicroActuators

The electrostatic actuator has lengtha, width b, and thickness h. The strain energy of the actuator is [25]
(4)U=12∫0b∫0aD(∂2w(x,y,t)∂x2+∂2w(x,y,t)∂y2)2dxdy,
where w is the deflection of the actuator, D=Eh3/(12(1-ν2)) is the flexural rigidity, E is Young's modulus, ν is the Poisson's ratio, and h is the actuator thickness. The kinetic energy of the microactuator is
(5)T=12ρh∫0b∫0a(∂w(x,y,t)∂t)2dxdy,
where t is the time and ρ is the density of the actuator material. Equations (4) and (5) are substituted into Hamilton equation as follows:
(6)∫t1t2(δT-δU+δW)dt=0,
where δW is the virtual work. This leads to the following equations for electrostatic actuator motion:
(7)D∂4w(x,y,t)∂x4+2D∂4w(x,y,t)∂x2∂y2+D∂4w(x,y,t)∂y4+ρh∂2w(x,y,t)∂t2=0,
The boundary conditions of the actuator are as follows:
(8)w(0,0,t)=0,w(a,0,t)=0,w(0,b,t)=0,w(a,b,t)=0,∂2w(x,y,t)∂x2+ν∂2w(x,y,t)∂y2=0,forx=0,0<y<b,∂3w(x,y,t)∂x3+(2-ν)∂3w(x,y,t)∂x2∂y=0,forx=0,0<y<b,∂2w(x,y,t)∂x2+ν∂2w(x,y,t)∂y2=0,forx=a,0<y<b,∂3w(x,y,t)∂x3+(2-ν)∂3w(x,y,t)∂x2∂y=0,forx=a,0<y<b,∂2w(x,y,t)∂y2+ν∂2w(x,y,t)∂x2=0,for0<x<a,y=0,∂3w(x,y,t)∂y3+(2-ν)∂3w(x,y,t)∂x∂y2=0,for0<x<a,y=0,∂2w(x,y,t)∂y2+ν∂2w(x,y,t)∂x2=0,for0<x<a,y=b,∂3w(x,y,t)∂y3+(2-ν)∂3w(x,y,t)∂x∂y2=0,for0<x<a,y=b.
After substituting w(x,y,t)=W(x,y)eiωt into (7), (7) can be rewritten as follows:
(9)D∂4W(x,y)∂x4+2D∂4W(x,y)∂x2∂y2+D∂4W(x,y)∂y4=ω2ρhW(x,y),
where ω is the natural frequency of the actuator. The boundary conditions of the electrostatic actuator are as follows:
(10)W(0,0)=0,W(a,0)=0,W(0,b)=0,W(a,b)=0,∂2W(x,y)∂x2+ν∂2W(x,y)∂x2=0,forx=0,0<y<b,∂3W(x,y)∂x3+(2-ν)∂3W(x,y)∂x2∂y=0,forx=0,0<y<b,∂2W(x,y)∂x2+ν∂2W(x,y)∂y2=0,forx=a,0<y<b,∂3W(x,y)∂x3+(2-ν)∂3W(x,y)∂x2∂y=0,forx=a,0<y<b,∂2W(x,y)∂y2+ν∂2W(x,y)∂x2=0,for0<x<a,y=0,∂3W(x,y)∂y3+(2-ν)∂3W(x,y)∂x∂y2=0,for0<x<a,y=0,∂2W(x,y)∂y2+ν∂2W(x,y)∂x2=0,for0<x<a,y=b,∂3W(x,y)∂y3+(2-ν)∂3W(x,y)∂x∂y2=0,for0<x<a,y=b.


## 4. Finite-element Method Model

The commercially available FEMLAB software package is used to evaluate dynamic problems based on partial differential equations. To derive finite formulations, the following virtual work principle must be utilized in the following equations [23, 24]:
(11)δΠ=∫V(δ[ε]T[σ]-δ[w*]T[f])dV=0,
where [ε] is a strain matrix, [σ] is a stress matrix, [w*] is a displacement matrix, and [f] is a inertia force matrix. The finite-element method can convert a differential equation into a set of algebraic equations, assume the shape of the solution in the element domain, and satisfy equilibrium. The stress-strain relationship for linear conditions is as follows:
(12)[σ]=[D̅][ε],
where [D̅] is an elastic matrix. The variables are approximated with functions in the chosen finite-element spaces. The following finite-element solution is assumed [23, 24]:
(13)w*=[N̅][W*]cos⁡⁡(ωt),
where [N̅] is the matrix of any suitable assumed shape function and [W*] is the displacement matrix. The finite element assembles all elements to form a complete structure to equilibrate a structure with its environment. The equation for the finite-element model of the electrostatic microactuator is as follows:
(14)[K*][W*]=ω2[M*][W*],
where [M*] is the mass matrix and [K*] is the stiffness matrix. The stiffness matrix can be written as follows:
(15)[K*]=∫V[[L̅][N̅]]T[D̅][[L̅][N̅]]dV,
where [L̅] is a linear differential operator matrix. The mass matrix can be written as follows:
(16)[M*]=∫V[N̅]T[ρ][N̅]dV.
The matrix order is N*×N*, where N* is the number of nodes for which the solution is unknown. Assembly of all element stiffness matrices and element mass matrices of the electrostatic actuator defines the following eigenvalue equation:
(17)|[K*]-ω2[M*]|=0.
The eigenvalues of the electrostatic actuator can be derived from (17) and are known as eigenfrequency problems.


## 5. Differential Quadrature Formula

The vibration response of the microactuator is numerically modeled using the differential quadrature method in this work. The differential quadrature method is used to convert the partial differential equations of the plates into a discrete eigenvalue problem. The roots of shifted Chebyshev and Legendre sampling point equation are used to select the sampling points in these analyses. The integrity and computational efficiency of the differential quadrature method in this problem is demonstrated below in several case studies. The differential quadrature method is a relatively new method that was introduced by Bellman and Casti [26]. After its appearance, several researchers have applied the differential quadrature method to solve a variety of problems in different fields of science and engineering. The differential quadrature method has been shown to be a powerful contender in solving initial and boundary value problems. Bert et al. [27–30] solved static and free vibration analysis of beams, plates, and compressible lubrication using the differential quadrature method. Chen and Zhong [31] reported that, due to their global domain properties, differential quadrature and differential cubature methods could solve nonlinear problems more efficiently than traditional numerical techniques such as the finite-element and the finite-difference methods. Civan and Sliepcevich [32] solved multivariable mathematical models using the quadrature method and the cubature method. Han and Liew [33] analyzed the axisymmetric free vibration of moderately thick annular plates using the differential quadrature methodology. Xu and Mazumder [34] derived the rational ABCD matrix representing the high-speed interconnect using the differential quadrature method. The differential quadrature method assumes that the derivative of a function at a sampling point can be approximated as a weighted linear combination of the functional values at all of the sampling points in the domain. The number of equations is dependent upon the selected number of the sampling points. A differential quadrature approximation at the ith discrete point on a grid in the x-axis direction can be approximated by [29]
(18)∂mW(xi,y)∂xm=∑j=1NxA¯ij(m)W(xj,y),fori=1,2,…,Nx.
A differential quadrature approximation at the ith discrete point on a grid in the direction of y-axis direction may be written as
(19)∂mW(x,yi)∂ym=∑j=1NyB¯ij(m)W(x,yj),fori=1,2,…,Ny,
where A¯ij(m) and B¯ij(m) are the differential weighting coefficients. The test function can be written as
(20)W(x,y)=α(x)β(y),(21)α(x)=M(x)(x-xk)M*(xk),fork=1,2,…,Nx,(22)M(x)=∏m=1Nx(x-xm),(23)M*(xk)=∏m=1m≠kNx(xk-xm),fork=1,2,…,Nx,(24)β(x)=L(y)(y-yl)L*(yl),forl=1,2,…,Ny,(25)L(y)=∏m=1Ny(y-ym),(26)L*(yl)=∏m=1m≠lNy(yl-ym),forl=1,2,…,Ny.
After substituting (20) to (18) and (19), differential weighting coefficients are given as follows:(27)A¯im(1)=M*(xi)(xi-xm)M*(xm),fori≠m,i=1,2,…,Nx,m=1,2,…,Nx,A¯ii(1)=-∑m=1m≠iNxA¯im(1),B¯im(1)=L*(yi)(yi-ym)L*(xm),fori≠m,i=1,2,…,Ny,m=1,2,…,Ny,B¯ii(1)=-∑m=1m≠iNyB¯im(1).
The weighting matrix A(1) is [27–30]
(28)A(1)=[[A¯im(1)][0]⋯[0][0][A¯im(1)]⋯⋮⋮⋮⋱[0][0][0]⋯[A¯im(1)]]NxNy×NxNy,fori=1,2,…,Nx,m=1,2,…,Nx.
The weighting matrix B(1) is [27–30]
(29)B(1)=[B̅11(1)[I]B̅12(1)[I]⋯B̅1Ny(1)[I]B̅21(1)[I]B̅22(1)[I]⋯B̅2Ny(1)[I]⋮⋮⋱⋮B̅Ny1(1)[I]B̅Ny2(1)[I]⋯B̅NyNy(1)[I]]NxNy×NxNy,
where [I] is an identity matrix of dimension Nx. The A(1) and B(1) are both square matrices of dimension NxNy. The higher-order derivates may be obtained using the following equations:
(30)Ai,j(m)=∑k=1NxNyAi,k(1)Ak,j(m-1),Bi,j(n)=∑k=1NxNyBi,k(1)Bk,j(n-1),A(m)Bi,j(n)=∑k=1NxNyAi,k(m)Bk,j(n).
The above relation gives the higher order weighting coefficient matrix based on the first-order derivative weighting coefficients. The selection of locations of the sampling points is important for ensuring the accuracy of the solution of differential equations. Using equally spaced points can be considered a convenient and easy selection method. The domain is divided by Nx×Ny points. The equally spaced sampling points are [29]
(31)xi=i-1Nx-1,fori=1,2,…,Nx,
for x direction, and(32)yi=i-1Ny-1,fori=1,2,…,Ny,
for y direction. An accurate solution can be obtained by choosing a set of unequally spaced sampling points. A simple and effective choice is the roots of shifted Chebyshev and Legendre points. Bert et al. [27–30] demonstrated that the following nonuniform grid spacing gives better and more reliable calculation results. The inner points are
(33)xi=12[1-cos⁡⁡(i-1)πNx-1],fori=1,2,…,Nx,
for x direction,
(34)yi=12[1-cos⁡⁡(i-2)πNy-3],fori=1,2,…,Ny,
for y direction. The use of zeros of shifted Legendre polynomials have been known to give good results. Although a convenient and commonly used choice of sampling of a quadrature grid is equally spaced points, nonuniformly spaced points generally achieve better accuracy than a quadrature solution does. Equations (18) and (19) are substituted for (9) and (10) as follows: (35)[DAi,1(4)a4+2DA(2)Bi,1(2)a2b2+DBi,1(4)b4DAi,2(4)a4+2DA(2)Bi,2(2)a2b2+DBi,2(4)b4⋯DAi,NxNy(4)a4+2DA(2)Bi,NxNy(2)a2b2+DBi,NxNy(4)b4]×[W1W2⋯WNxNy]T=[ω2ρhWi],fori=1,2,3,…,NxNy-1,NxNy,[Wi]=[0],fori=1,Nx,(Ny-1)Nx+1,NxNy,[Ai,1(2)a2+vBi,1(2)b2Ai,2(2)a2+vBi,2(2)b2⋯Ai,NxNy(2)a2+vBi,NxNy(2)b2][W1W2⋯WNxNy]T=[0],fori=Nx+1,2Nx+1,3Nx+1,4Nx+1,…,(Ny-3)Nx+1,(Ny-2)Nx+1,2Nx,3Nx,4Nx,5Nx,…,(Ny-2)Nx,(Ny-1)Nx,[Ai,1(3)a3+(2-v)A(2)Bi,1(1)a2bAi,2(3)a3+(2-v)A(2)Bi,2(1)a2b⋯Ai,NxNy(3)a3+(2-v)A(2)Bi,NxNy(1)a2b]×[W1W2⋯WNxNy]T=[0],fori=Nx+1,2Nx+1,3Nx+1,4Nx+1,…,(Ny-3)Nx+1,(Ny-2)Nx+1,2Nx,3Nx,4Nx,5Nx,…,(Ny-2)Nx,(Ny-1)Nx,[Bi,1(2)b2+vAi,1(2)a2Bi,2(2)b2+vAi,2(2)a2⋯Bi,NxNy(2)b2+vAi,NxNy(2)a2][W1W2⋯WNxNy]T=[0],fori=2,3,4,5,…,Nx-2,Nx-1,(Ny-1)Nx+2,(Ny-1)Nx+3,(Ny-1)Nx+4,(Ny-1)Nx+5,…,NxNy-2,NxNy-1,[Bi,1(3)b3+(2-v)B(2)Ai,1(1)b2aBi,2(3)b3+(2-v)B(2)Ai,2(1)b2a⋯Bi,NxNy(3)b3+(2-v)B(2)Ai,NxNy(1)b2a]×[W1W2⋯WNxNy]T=[0],fori=2,3,4,5,…,Nx-2,Nx-1,(Ny-1)Nx+2,(Ny-1)Nx+3,(Ny-1)Nx+4,(Ny-1)Nx+5,…,NxNy-2,NxNy-1.Equation (35) is solved to obtain the frequencies of the electrostatic microactuators.


## 6. Frequencies of Electrostatic Micro Actuators

The material parameters of the electrostatic actuator areρ=2.328×103 kg/m3 and E=150 GPa [35]. The FEMLAB finite-element package is used to analyze the model. The 26733 solid elements are used to calculate the frequency of the actuator. Figure 6 plots the frequencies of the microelectrostatic actuator with various lengths. The width and thickness of the electrostatic actuator are 176 and 8 micrometers, respectively. The numerical results indicate that the frequencies calculated using the finite-element method and the differential quadrature method are almost identical. Several case studies have validated the applicability of the method for solving such engineering problems. The method provides accurate results with relatively less computational and modeling efforts. The calculation results show that higher lengths produce smaller frequencies of the electrostatic actuator. The differential quadrature method has become a preferred method to the finite-element method. Figure 7 plots the frequencies of the electrostatic actuator with various widths. The length and thickness of the electrostatic actuator are 262 and 8 micrometers, respectively. The numerical results in this example show that the widths can significantly affect the dynamic behavior of the electrostatic actuator. Higher widths produce smaller frequencies of the electrostatic actuator. Figure 8 shows the frequencies of the microelectrostatic actuator with various thicknesses. The length and width of the electrostatic actuator are 262 and 176 micrometers, respectively. The numerical results indicated that the frequency of the microactuator is increased for the actuator with a larger value of thickness. The numerical results indicate that the thickness of the actuator is a very sensitive parameter to the frequency of the actuator. The frequency of the actuator increases with increases in thickness h.Figure 6
The lowest six frequencies of electrostatic actuator for various lengths.Figure 7
The lowest six frequencies of electrostatic actuator for various widths.Figure 8
The lowest six frequencies of electrostatic actuator for various thicknesses.


## 7. Conclusions

Numerical results indicate that length and width significantly impact the capacity of electrostatic microactuator. The presented formulation reveals that the differential quadrature approach is convenient for solving problems governed by fourth- or higher-order differential equations. Simulation results verify that the differential quadrature method obtains accurate results with relatively minimal computational and modeling efforts. Length, width, and thickness can markedly affect electrostatic microactuator frequencies. The FEMLAB can handle capacitance and dynamic problems as well. The differential quadrature methodology may be further examined to solve more complicated problems or in other fields of science.



---

*Source: 576079-2011-08-02.xml*