# A Wavelet Interpolation Galerkin Method for the Simulation of MEMS Devices under the Effect of Squeeze Film Damping

Authors: Pu Li; Yuming Fang
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

This paper presents a new wavelet interpolation Galerkin method for the numerical simulation of MEMS devices under the effect of squeeze film damping. Both trial and weight functions are a class of interpolating functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions. To the best of our knowledge, this is the first time that wavelets have been used as basis functions for solving the PDEs of MEMS devices. As opposed to the previous wavelet-based methods that are all limited in one energy domain, the MEMS devices in the paper involve two coupled energy domains. Two typical electrically actuated micro devices with squeeze film damping effect are examined respectively to illustrate the new wavelet interpolation Galerkin method. Simulation results show that the results of the wavelet interpolation Galerkin method match the experimental data better than that of the finite difference method by about 10%.


---

## Body

## 1. Introduction

Modeling and simulation of MEMS devices play an important role in the design phase for system optimization and for the reduction of design cycles. The performances of MEMS devices are represented by partial-differential equations (PDEs) and associated boundary conditions. In the past two decades, there have been extensive, and successful, works focused on solving the partial-differential equations of MEMS [1–15]. A detailed review of the works is available in [1]. In the previous works, Galerkin method was widely used to reduce the partial-differential equations to ordinary-differential equations (ODEs) in time and then solve the reduced equations either numerically or analytically. The previous works differ from each other in the choice of the basis functions.The basis set can be chosen arbitrarily, as long as its elements satisfy all of the boundary conditions and are sufficiently differentiable. To enhance convergence, the basis set has to be chosen to resemble the behavior of the device. For example, two ways have been used to generate the basis set for the reduced-order models of MEMS devices [1]. The first way [4, 9] uses the undamped linear model shapes of the undeflected microstructure as basis functions. For simple structures with simple boundary conditions, the mode shapes are found analytically. For complex structures or complex boundary conditions, the linear mode shapes are obtained numerically using the finite element method. The second way [2] conducts experiments or solves the PDEs using FEM or FDM to generate snapshots under a training signal, then applies a modal analysis method (one of the variation of the proper orthogonal decomposition method [6]) to the time series to extract the mode shapes of the device structural elements.In the past two decades also, a new numerical concept was introduced and is gaining increasing popularity [16–25]. The method is based on the expansion of functions in terms of a set of basis functions called wavelets. Indeed wavelets have many excellent properties such as orthogonality, compact support, exact representation of polynomials to a certain degree, and flexibility to represent functions at different levels of resolution. Indeed a complete basis can be generated easily by a signal function through dilatation and translation. The wavelet-based methods may be classified as wavelet-Galerkin method [19, 20], wavelet-collocation method [21, 22], and wavelet interpolation Galerkin method [23–25]. Among the three methods, the wavelet-Galerkin method is the most common one because of its implementation simplicity. The method is a Galerkin scheme using scaling or wavelet functions as the trial and weight functions. However, both scaling and wavelet functions do not satisfy the boundary conditions. Thus the treatment of general boundary conditions is a major difficulty for the application of the wavelet-Galerkin method, especially for the bounded region problems, even though different efforts [19, 20] have been made. For the wavelet-collocation method, boundary conditions can be treated in a satisfactory way [21]. In the method, trial functions are a class of interpolating functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions. However, the method requires the calculation of higher-order derivatives (up to the second derivatives for second-order parabolic problems) of the wavelets. Due to the derivatives of compactly supported wavelets being highly oscillatory, it is difficult to compute the connection coefficients by the numerical evaluation of integral [18]. The wavelet interpolation Galerkin method is a Galerkin scheme that both trial and weight functions are a class of interpolating functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions. For the method, the boundary conditions [24] can be treated easily and the formulations are derived from the weak form; thus only the first derivatives of wavelets (for second-order parabolic problems) are required.Wavelets have proven to be an efficient tool of analysis in many fields including the solution of PDEs. However, few papers in MEMS area give attention to the wavelet-based methods. This paper presents a new wavelet interpolation Galerkin method for the numerical simulation of MEMS devices under the effect of squeeze film damping. To the best of our knowledge, this is the first time that wavelets have been used as basis functions for solving the PDEs of MEMS devices. As opposed to the previous wavelet-based methods that are all limited in one energy domain, the MEMS devices in the paper involve two coupled energy domains. The squeeze film damping effect on the dynamics of microstructures has already been extensively studied. We stress that our intention here is not to discover new physics to the squeeze film damping.The outline of this paper is as follows. Section2 presents a brief introduction to some major concepts and properties of wavelets. In Sections 3 and 4, two typical electrically actuated micro devices with squeeze film damping effect are examined respectively to illustrate the wavelet interpolation Galerkin method. Section 5 calculates the frequency responses and the quality factors using the present method, and compares the calculated results with those generated by experiment [26, 27], by the finite difference method, and by other published analytical models [15, 26]. Finally, a conclusion is given in Section 6.


## 2. Basic Concepts of Daubechies’ Wavelets and Wavelet Interpolation

In this section, we shall give a brief introduction to the concepts and properties of Daubechies’ wavelets. More detailed discussions can be found in [16–18, 21].


### 2.1. Daubechies’ Orthonormal Wavelets

Daubechies [16, 17] constructed a family of orthonomal bases of compactly supported wavelets for the space of square-integrable funcntions, L2(R). Due to the fact that they possess several useful properties, such as orthogonality, compact support, exact representation of polynomials to a certain degree, and ability to represent functions at different levels of resolution, Daubechies’ wavelets have gained great interest in the numerical solutions of PDEs [18–22].Daubechies’ functions are easy to construct [16, 17]. For an even integer L, we have the Daubechies’ scaling function ϕ(x) and wavelet ψ(x) satisfying(2.1)ϕ(x)=∑i=0L-1p̂iϕ(2x-i)ψ(x)=∑i=2-L1(-1)ip̂1-iϕ(2x-i).
The fundamental support of the scaling function ϕ(x) is in the interval [0,L-1] while that of the corresponding wavelet ψ(x) is in the interval [1-L/2,L/2]. The parameter L will be referred to as the degree of the scaling function ϕ(x). The coefficients p̂i are called the wavelet filter coefficients. Daubenchies  [16, 17] established these wavelet filter coefficients to satisfy the following conditions:(2.2)∑i=0L-1p̂i=2,∑i=0L-1p̂ip̂i-m=δ0,m,∑i=2-L1(-1)ip̂1-ip̂i-2m=0forinteger,∑i=0L-1(-1)iimp̂i=0,m=0,1,…,L2-1,
where δ0,m is the Kronecker delta function. Correspondingly, the constructed scanling function ϕ(x) and wavelet ψ(x) have the following properties:(2.3)∫-∞∞ϕ(x)dx=1,∫-∞∞ϕ(x-i)ϕ(x-m)dx=δi,mforintegersi,m,∫-∞∞ϕ(x)ψ(x-m)dx=0forintegerm,∫-∞∞xmψ(x)dx=0,m=0,1,…,L2-1.Denote byL2(R) the space of square-integrable functions on the real line. Let VJ and WJ be the subspace generated, respectively, as the L2-closure of the linear spans of ϕJ,i(x)=2J/2ϕ(2Jx-i) and ψJ,i(x)=2J/2ψ(2Jx-i), J,i∈Z. ℤ denotes the set of integers. Then (2.3) implies that
(2.4)VJ+1=VJ⊕WJ,V0⊂V1⊂⋯VJ⊂VJ+1,VJ+1=V0⊕W0⊕W1⊕⋯WJ,
Equation (2.4) presents the multiresolution properties of    wavlets. Any function f∈L2(R), may be approximated by the multiresolution apparatus described above, by its projection PJVf onto the subspaceVJ(2.5)PJVf=∑i∈ZfJ,iϕJ,i(x).


### 2.2. Wavelet Interpolation Scaling Function

For a given Daubechies’ scaling function, its autocorrelation functionθ(x) can be defined as follows [21]:
(2.6)θ(x)=∫-∞∞ϕ(τ)ϕ(τ-x)dτ.
The function satisfies the following interpolating property(2.7)θ(k)=δ0,k,k∈Z,
and has a symmetric support [-(L-1),(L-1)]. The derivative of the function θ(k) may be computed by differentiating the convolution product(2.8)θ(s)(k)=(-1)(s)∫-∞∞ϕ(τ)ϕ(s)(τ-k)dτ.
Let θ(x) act as the scaling function, we have(2.9)θJ,k(x)=θ(2Jx-k),k∈Z.
For a set of dyadic grids of the type xkJ∈R: xkJ=2-Jk, where k,J∈Z, the θJ,k(x) verifies the interpolation property at the dyadic points: θJ,k(xnJ)=δn,k. Let VJx be the linear span of the set {θ(2Jx-k),k∈Z}. It can be proved that {VJx} forms a multiresolution analysis, where θJ,k(x) acts as the role of scaling function (the so-called interpolation scaling function), and the set {θ(2Jx-k),k∈Z} is a Riesz’s basis for VJx. For a function f∈H1(R), an interpolation operator IJ:H1(R)→VJx can be defined [21]:
(2.10)IJ(f)=∑kfkJθ(2Jx-k),k∈Z,
where fkJ=f(xkJ)=f(2-Jk). Thus, for a function f(x) defined on x∈[0,1], f(x) has the following approximation(2.11)f(x)=∑k=-(L-1)2J+(L-1)fJ,kθ(2Jx-k)=∑k=-(L-1)1fJ,kθ(2Jx-k)+∑k=02JfJ,kθ(2Jx-k)+∑k=2J+12J+L-1fJ,kθ(2Jx-k).In this paper, wavelet collocation scheme is applied onx∈[0,1], where xkJ=2-Jk and k=0,1,…,2J. Therefore, instead of the values of f(x) at xkJ, k=-(L-1),…,-1 and k=2J+1,…,(L-1), we may use some values which are extrapolated from the values in those dyadic points internal to the interval x∈[0,1]. As described in [21, 22], we define(2.12)f(x)=∑k=02JfJ,kθ̂(2Jx-k),
where(2.13)θ̂(2Jx-k)={θ(2Jx-k)+∑n=-(L-1)-1ankθ(2Jx-n),k=0,1,…,2M-1θ(2Jx-k),k=2M,…,2J-2Mθ(2Jx-k)+∑n=2J+12J+L-1bnkθ(2Jx-n),k=2J-2M+1,…,2J,
where the coefficients ank and bnk are defined by(2.14)ank=lk1(xnJ),bnk=lk2(xnJ),
where lk1(x) and lk2(x) represent Lagrange interpolation polynomials, defined by(2.15)lk1(x)=∏i=0i≠k2M-1x-xiJxkJ-xiJ,lk2(x)=∏i=2J-2M+1i≠k2Jx-xiJxkJ-xiJ.An analogous manner can be given for two-dimensional problem. By using tensor products, it is then possible to define a multiresolution on the squarex,y∈[0,1]. The two-dimensional scaling function is defined by Θk,k'J(x,y)=∑k=02J∑k'=02Jfk,k'Jθ̂(2Jx-k)θ̂(2Jy-k'). Let VJxy=VJx⊗VJy be the linear span of the set {θ(2Jx-k)θ(2Jy-k′),J,k,k'∈Z}; thus the set {VJxy} forms a multiresolution analysis and the set {θ(2Jx-k)θ(2Jy-k′),k,k'∈Z} is a Riesz basis for {VJxy}. Therefore, for a function f(x,y) defined on x,y∈[0,1], it has the following approximation:(2.16)f(x,y)=∑k=02J∑k'=02Jfk,k'Jθ̂(2Jx-k)θ̂(2Jy-k').


## 2.1. Daubechies’ Orthonormal Wavelets

Daubechies [16, 17] constructed a family of orthonomal bases of compactly supported wavelets for the space of square-integrable funcntions, L2(R). Due to the fact that they possess several useful properties, such as orthogonality, compact support, exact representation of polynomials to a certain degree, and ability to represent functions at different levels of resolution, Daubechies’ wavelets have gained great interest in the numerical solutions of PDEs [18–22].Daubechies’ functions are easy to construct [16, 17]. For an even integer L, we have the Daubechies’ scaling function ϕ(x) and wavelet ψ(x) satisfying(2.1)ϕ(x)=∑i=0L-1p̂iϕ(2x-i)ψ(x)=∑i=2-L1(-1)ip̂1-iϕ(2x-i).
The fundamental support of the scaling function ϕ(x) is in the interval [0,L-1] while that of the corresponding wavelet ψ(x) is in the interval [1-L/2,L/2]. The parameter L will be referred to as the degree of the scaling function ϕ(x). The coefficients p̂i are called the wavelet filter coefficients. Daubenchies  [16, 17] established these wavelet filter coefficients to satisfy the following conditions:(2.2)∑i=0L-1p̂i=2,∑i=0L-1p̂ip̂i-m=δ0,m,∑i=2-L1(-1)ip̂1-ip̂i-2m=0forinteger,∑i=0L-1(-1)iimp̂i=0,m=0,1,…,L2-1,
where δ0,m is the Kronecker delta function. Correspondingly, the constructed scanling function ϕ(x) and wavelet ψ(x) have the following properties:(2.3)∫-∞∞ϕ(x)dx=1,∫-∞∞ϕ(x-i)ϕ(x-m)dx=δi,mforintegersi,m,∫-∞∞ϕ(x)ψ(x-m)dx=0forintegerm,∫-∞∞xmψ(x)dx=0,m=0,1,…,L2-1.Denote byL2(R) the space of square-integrable functions on the real line. Let VJ and WJ be the subspace generated, respectively, as the L2-closure of the linear spans of ϕJ,i(x)=2J/2ϕ(2Jx-i) and ψJ,i(x)=2J/2ψ(2Jx-i), J,i∈Z. ℤ denotes the set of integers. Then (2.3) implies that
(2.4)VJ+1=VJ⊕WJ,V0⊂V1⊂⋯VJ⊂VJ+1,VJ+1=V0⊕W0⊕W1⊕⋯WJ,
Equation (2.4) presents the multiresolution properties of    wavlets. Any function f∈L2(R), may be approximated by the multiresolution apparatus described above, by its projection PJVf onto the subspaceVJ(2.5)PJVf=∑i∈ZfJ,iϕJ,i(x).


## 2.2. Wavelet Interpolation Scaling Function

For a given Daubechies’ scaling function, its autocorrelation functionθ(x) can be defined as follows [21]:
(2.6)θ(x)=∫-∞∞ϕ(τ)ϕ(τ-x)dτ.
The function satisfies the following interpolating property(2.7)θ(k)=δ0,k,k∈Z,
and has a symmetric support [-(L-1),(L-1)]. The derivative of the function θ(k) may be computed by differentiating the convolution product(2.8)θ(s)(k)=(-1)(s)∫-∞∞ϕ(τ)ϕ(s)(τ-k)dτ.
Let θ(x) act as the scaling function, we have(2.9)θJ,k(x)=θ(2Jx-k),k∈Z.
For a set of dyadic grids of the type xkJ∈R: xkJ=2-Jk, where k,J∈Z, the θJ,k(x) verifies the interpolation property at the dyadic points: θJ,k(xnJ)=δn,k. Let VJx be the linear span of the set {θ(2Jx-k),k∈Z}. It can be proved that {VJx} forms a multiresolution analysis, where θJ,k(x) acts as the role of scaling function (the so-called interpolation scaling function), and the set {θ(2Jx-k),k∈Z} is a Riesz’s basis for VJx. For a function f∈H1(R), an interpolation operator IJ:H1(R)→VJx can be defined [21]:
(2.10)IJ(f)=∑kfkJθ(2Jx-k),k∈Z,
where fkJ=f(xkJ)=f(2-Jk). Thus, for a function f(x) defined on x∈[0,1], f(x) has the following approximation(2.11)f(x)=∑k=-(L-1)2J+(L-1)fJ,kθ(2Jx-k)=∑k=-(L-1)1fJ,kθ(2Jx-k)+∑k=02JfJ,kθ(2Jx-k)+∑k=2J+12J+L-1fJ,kθ(2Jx-k).In this paper, wavelet collocation scheme is applied onx∈[0,1], where xkJ=2-Jk and k=0,1,…,2J. Therefore, instead of the values of f(x) at xkJ, k=-(L-1),…,-1 and k=2J+1,…,(L-1), we may use some values which are extrapolated from the values in those dyadic points internal to the interval x∈[0,1]. As described in [21, 22], we define(2.12)f(x)=∑k=02JfJ,kθ̂(2Jx-k),
where(2.13)θ̂(2Jx-k)={θ(2Jx-k)+∑n=-(L-1)-1ankθ(2Jx-n),k=0,1,…,2M-1θ(2Jx-k),k=2M,…,2J-2Mθ(2Jx-k)+∑n=2J+12J+L-1bnkθ(2Jx-n),k=2J-2M+1,…,2J,
where the coefficients ank and bnk are defined by(2.14)ank=lk1(xnJ),bnk=lk2(xnJ),
where lk1(x) and lk2(x) represent Lagrange interpolation polynomials, defined by(2.15)lk1(x)=∏i=0i≠k2M-1x-xiJxkJ-xiJ,lk2(x)=∏i=2J-2M+1i≠k2Jx-xiJxkJ-xiJ.An analogous manner can be given for two-dimensional problem. By using tensor products, it is then possible to define a multiresolution on the squarex,y∈[0,1]. The two-dimensional scaling function is defined by Θk,k'J(x,y)=∑k=02J∑k'=02Jfk,k'Jθ̂(2Jx-k)θ̂(2Jy-k'). Let VJxy=VJx⊗VJy be the linear span of the set {θ(2Jx-k)θ(2Jy-k′),J,k,k'∈Z}; thus the set {VJxy} forms a multiresolution analysis and the set {θ(2Jx-k)θ(2Jy-k′),k,k'∈Z} is a Riesz basis for {VJxy}. Therefore, for a function f(x,y) defined on x,y∈[0,1], it has the following approximation:(2.16)f(x,y)=∑k=02J∑k'=02Jfk,k'Jθ̂(2Jx-k)θ̂(2Jy-k').


## 3. Wavelet Interpolation Galerkin Method for a Parallel Plate   Microresonator under the Effect of Squeeze Film Damping

### 3.1. Governing Equations

In this section, we examine the example of a rectangular parallel plate under the effect of squeeze film damping. As shown in Figure1, the rectangular parallel plate is excited by a conventional voltage. The voltage is composed of a dc component V0 and a small ac component v(t), V0≫v(t). The plate is rigid. The displacement of the plate under the electric force is composed of a static component to the dc voltage, denoted by z0, and a small dynamic component due to the ac voltage, denoted by z(t), z0≫z(t), that is,
(3.1)zE(t)=z0+z(t).
The equation of motion that governs the displacement of the plate is written as(3.2)mplatez̈E+kspringzE=εAplate(V0+v)22(g0-zE)2-f(t),
where mplate is the mass of the plate, Aplate is the are of the plate, kspring is the stiffness of the spring, g0 is the zero-voltage air gap spacing, ε is the dielectric constant of the gap medium, f(t) is the force acting on the plate owing to the pressure of the squeeze gas film between the plate and the substrate.A schematic drawing of an electrically actuated microplate under the effect of squeeze film damping.
(a)
Side view of the microplate(b)
Top view of the microplateWe expand (3.2) in a Taylor series around V0 and z0 up to first order and rewrite (3.2) as(3.3)mplatez̈+kEz=εAplateV0ĝ02v-f(t),
where kE=kspring-(εAplateV02/(g0-z0)3), ĝ0=g0-z0. The force f(t) acting on the plate owing to the pressure of the squeeze gas film is given by
(3.4)f(t)=∫0ly∫0lx(p(x,y,t)-p0)dxdy,
where lx and ly are the length and width of plate, p(x,y,t) is the absolute pressure in the gap and p0 is the ambient pressure. The pressure p(x,y,t) is governed by the nonlinear Reynolds equation [3](3.5)∂∂x(h3p∂p∂x)+∂∂y(h3p∂p∂y)=12ηeff(h∂p∂t+p∂h∂t),
where h(x,t)=g0-z0-z(t)=ĝ0-z(t) and ηeff is the effective viscosity of the fluid in the gap. In this section, all edges of the rectangular plate are ideally vented; thus the pressure boundary conditions for the case in Figure 1 are(3.6)p(x,0,t)=p(x,ly,t)=p(0,y,t)=p(lx,y,t)=p0.For convenience, we introduce the nondimensional variables(3.7)X=xlx,Y=yly,Z=zĝ0,P̂=pp0,T=tS,H=hĝ0=1-Z,
where T is a timescale, S=mplate/kE=1/ωn, ωn is the nature frequency of the plate. Substituting (3.7) into (3.3)–(3.6), we obtain(3.8)d2ZdT2+Z=αV0v-Pnon∬01(P̂-1)dXdY,(3.9)∂∂X(H3P̂∂P̂∂X)+β2∂∂Y(H3P̂∂P̂∂Y)=σS(H∂P̂∂T+P̂∂H∂T),
where α=εAplate/kEĝ03, Pnon=p0lxly/kEĝ0, σ=12ηefflx2/ĝ02p0, and β=lx/ly. The nondimensional boundary conditions are(3.10)P̂(X,0,T)=P̂(X,1,T)=P̂(1,Y,T)=P̂(0,Y,T)=1.As mentioned above, the microplate is under small oscillation aroundĝ0 and therefore the pressure variation from ambient in the squeeze film is also small, P̂(X,Y,T) is given by(3.11)P̂(X,Y,T)=pp0=1+P(X,Y,T),
where |P(X,Y,T)|≪1. Substituting (3.11) into (3.9), and linearizing the outcome around p0 and ĝ0, we obtain
(3.12)∂2P∂X2+β2∂2P∂Y2-σS∂P∂T=-σS∂Z∂T.
The boundary conditions for the case are(3.13)P(X,0,T)=P(X,1,T)=P(0,Y,T)=P(1,Y,T)=0.For a harmonic excitation, the ac component voltagev(t) is given by(3.14)v(T)=v0ejωTS.
Usually, the excitation frequency ω is approximate to the natural frequency ωn. The steady-state solution of (3.8) and (3.12) may be expressed by(3.15)Z(T)=AejωTS,(3.16)P(X,Y,T)=A·PA(X,Y)ejωTS,
where A is the complex amplitude to be determined. Substituting (3.15) and (3.16) into (3.12), we obtain(3.17)∂2PA(X,Y)∂X2+β2∂2PA(X,Y)∂Y2-jσωPA(X,Y)=-jσω.
The boundary conditions are(3.18)PA(X,0)=PA(X,1)=PA(0,Y)=PA(1,Y)=0.


### 3.2. Wavelet Interpolation Method for Squeeze Film Damping Equations

#### 3.2.1. Construction of Basis Functions

In this subsection, the approximate solution ofPA(X,Y) is approximated by the following form:
(3.19)PA(X,Y)≈∑k=02J∑k′=02Jpk,k′JΘk,k′J(X,Y)=∑k=02J∑k′=02Jpk,k′Jθ̂J,k(X)θ̂J,k′(Y)=∑k=02J∑k′=02Jpk,k′Jθ̂(2JX-k)θ̂(2JY-k′),k,k′∈Z,
where the unknowns pk,k′J are the values of PA(X,Y) at the dyadic points X=k2-J, and  Y=k′2-J. The unknowns pk,k′J are complex.For the application of Galerkin method, (3.19) should be able to satisfy the boundary conditions. Substituting (3.19) into (3.18), leads to(3.20)∑k′=02Jp0,k′Jθ̂(0)θ̂(2JY-k′)=0⇒p0,k′J=0,fork′=0,1,2,…2J,∑k′=02Jp2J,k′Jθ̂(0)θ̂(2JY-k′)=0⇒p2J,k′J=0,fork′=0,1,2,…2J,∑k=12J-1pk,0Jθ̂(2JX-k)θ̂(0)=0⇒pk,0J=0,fork=1,2,…,(2J-1),∑k=12J-1pk,2JJθ̂(2JX-k)θ̂(0)=0⇒pk,2JJ=0,fork=1,2,…,(2J-1).
Thus (3.19) is rewritten as (3.21)PA(X,Y)=∑k=12J-1∑k′=12J-1pk,k′JΘk,k′J(X,Y)=∑k=12J-1∑k′=12J-1pk,k′Jθ̂(2JX-k)θ̂(2JY-k′).


#### 3.2.2. Discretion of the Boundary Value Problem

The weak form functional of (3.17) is(3.22)W(PA)=∬Ω{12[(∂PA∂X)2+β2(∂PA∂Y)2+jσωPA2]-jσωPA}dXdY.
From the necessary conditions for the determination of the minimum W, we obtain(3.23)δW(PA)=∬Ω[∂δPA∂X∂PA∂X+β2∂δPA∂Y∂PA∂Y+jσωδPAPA-jσωδPA]dXdY=0.
Substituting (3.19) into (3.23), leads to(3.24)∑k=12J-1∑k′=12J-1{∬Ω[∂Θm,nJ∂X∂Θk,k′J∂X+β2∂Θm,nJ∂Y∂Θk,k′J∂Y+jσωΘm,nJΘk,k′J]dXdY}pk,k′J=jσω∬ΩΘm,nJdXdY,form,n=1,2,…,(2J-1).
This is a (2J-1)2×(2J-1)2 linear system(3.25)Θp=jσωE,
where p=[p1,1Jp1,2J⋯p1,2J-1Jp2,1Jp2,2J⋯p2J-1,2J-1J]T is an (2J-1)2×1 unknown coefficients’ vector, E=[∬ΩΘ1,1JdXdY∬ΩΘ1,2JdXdY⋯∬ΩΘ2J-1,2J-1JdXdY]T is a (2J-1)2×1 matrix, and Θ is a (2J-1)2×(2J-1)2 matrix. The entries in Θ are of the form(3.26)Θ̂(i,j)=Θ̂([(m-1)(2J-1)+n],[(k-1)(2J-1)+k′])=∬Ω[∂θ̂(2JX-m)∂Xθ̂(2JY-n)∂θ̂(2JX-k)∂Xθ̂(2JY-k′)+β2θ̂(2JX-m)∂θ̂(2JY-n)∂Yθ̂(2JX-k)∂θ̂(2JY-k′)∂Y+jσωθ̂(2JX-m)θ̂(2JY-n)θ̂(2JX-k)θ̂(2JY-k′)∂θ̂(2JX-m)∂X]dXdY.


#### 3.2.3. Squeeze Film Damping of the Parallel Plate

The numerical solution of (3.25) can be written as(3.27)p=jσωΘ-1E.
The elements of p can be expressed as(3.28)pk,k′J=pk,k′J,R+jpk,k′J,Iform,n=1,2,…,(2J-1),
where pk,k′J,R and pk,k′J,I are the real and imaginary parts of pk,k′J, respectively. Using (3.21) and (3.28), the force acting on the plate owing to the pressure of the squeeze gas film can be rewritten as
(3.29)Pnon∬01(P̂-1)dXdY=AejωTS·Pnon∑k=12J-1∑k′=12J-1(pk,k′J,R+jpk,k′J,I)∬01θ̂(2JX-k)θ̂(2JY-k′)dXdY=Ka·Z(T)+Ca·dZ(T)dT,
where(3.30)Ka=Pnon∑k=12J-1∑k′=12J-1pk,k′J,R∬01θ̂(2JX-k)θ̂(2JY-k′)dXdY,Ca=PnonωS∑k=12J-1∑k′=12J-1pk,k′J,I∬01θ̂(2JX-k)θ̂(2JY-k′)dXdYKa·Z(T) and Ca·(dZ(T)/dT) are the spring and damping components of the force. Substituting (3.29), (3.14) and (3.15) into (3.8), we obtain(3.31)d2ZdT2+Ca·dZdT+(Ka+1)Z(T)=αV0v(T),Z(T)=AejωTS=aV0v0Ka+1·11-(ω2(S2/(Ka+1)))+(jω(CaS/(Ka+1)))ejωTS,
where S=mplate/kE=1/ωn. The quality factor and the damped natural frequency are expressed as(3.32)Qsqueeze=12ξ=Ka+1Ca,ωsqueeze=ωnKa+1.


## 3.1. Governing Equations

In this section, we examine the example of a rectangular parallel plate under the effect of squeeze film damping. As shown in Figure1, the rectangular parallel plate is excited by a conventional voltage. The voltage is composed of a dc component V0 and a small ac component v(t), V0≫v(t). The plate is rigid. The displacement of the plate under the electric force is composed of a static component to the dc voltage, denoted by z0, and a small dynamic component due to the ac voltage, denoted by z(t), z0≫z(t), that is,
(3.1)zE(t)=z0+z(t).
The equation of motion that governs the displacement of the plate is written as(3.2)mplatez̈E+kspringzE=εAplate(V0+v)22(g0-zE)2-f(t),
where mplate is the mass of the plate, Aplate is the are of the plate, kspring is the stiffness of the spring, g0 is the zero-voltage air gap spacing, ε is the dielectric constant of the gap medium, f(t) is the force acting on the plate owing to the pressure of the squeeze gas film between the plate and the substrate.A schematic drawing of an electrically actuated microplate under the effect of squeeze film damping.
(a)
Side view of the microplate(b)
Top view of the microplateWe expand (3.2) in a Taylor series around V0 and z0 up to first order and rewrite (3.2) as(3.3)mplatez̈+kEz=εAplateV0ĝ02v-f(t),
where kE=kspring-(εAplateV02/(g0-z0)3), ĝ0=g0-z0. The force f(t) acting on the plate owing to the pressure of the squeeze gas film is given by
(3.4)f(t)=∫0ly∫0lx(p(x,y,t)-p0)dxdy,
where lx and ly are the length and width of plate, p(x,y,t) is the absolute pressure in the gap and p0 is the ambient pressure. The pressure p(x,y,t) is governed by the nonlinear Reynolds equation [3](3.5)∂∂x(h3p∂p∂x)+∂∂y(h3p∂p∂y)=12ηeff(h∂p∂t+p∂h∂t),
where h(x,t)=g0-z0-z(t)=ĝ0-z(t) and ηeff is the effective viscosity of the fluid in the gap. In this section, all edges of the rectangular plate are ideally vented; thus the pressure boundary conditions for the case in Figure 1 are(3.6)p(x,0,t)=p(x,ly,t)=p(0,y,t)=p(lx,y,t)=p0.For convenience, we introduce the nondimensional variables(3.7)X=xlx,Y=yly,Z=zĝ0,P̂=pp0,T=tS,H=hĝ0=1-Z,
where T is a timescale, S=mplate/kE=1/ωn, ωn is the nature frequency of the plate. Substituting (3.7) into (3.3)–(3.6), we obtain(3.8)d2ZdT2+Z=αV0v-Pnon∬01(P̂-1)dXdY,(3.9)∂∂X(H3P̂∂P̂∂X)+β2∂∂Y(H3P̂∂P̂∂Y)=σS(H∂P̂∂T+P̂∂H∂T),
where α=εAplate/kEĝ03, Pnon=p0lxly/kEĝ0, σ=12ηefflx2/ĝ02p0, and β=lx/ly. The nondimensional boundary conditions are(3.10)P̂(X,0,T)=P̂(X,1,T)=P̂(1,Y,T)=P̂(0,Y,T)=1.As mentioned above, the microplate is under small oscillation aroundĝ0 and therefore the pressure variation from ambient in the squeeze film is also small, P̂(X,Y,T) is given by(3.11)P̂(X,Y,T)=pp0=1+P(X,Y,T),
where |P(X,Y,T)|≪1. Substituting (3.11) into (3.9), and linearizing the outcome around p0 and ĝ0, we obtain
(3.12)∂2P∂X2+β2∂2P∂Y2-σS∂P∂T=-σS∂Z∂T.
The boundary conditions for the case are(3.13)P(X,0,T)=P(X,1,T)=P(0,Y,T)=P(1,Y,T)=0.For a harmonic excitation, the ac component voltagev(t) is given by(3.14)v(T)=v0ejωTS.
Usually, the excitation frequency ω is approximate to the natural frequency ωn. The steady-state solution of (3.8) and (3.12) may be expressed by(3.15)Z(T)=AejωTS,(3.16)P(X,Y,T)=A·PA(X,Y)ejωTS,
where A is the complex amplitude to be determined. Substituting (3.15) and (3.16) into (3.12), we obtain(3.17)∂2PA(X,Y)∂X2+β2∂2PA(X,Y)∂Y2-jσωPA(X,Y)=-jσω.
The boundary conditions are(3.18)PA(X,0)=PA(X,1)=PA(0,Y)=PA(1,Y)=0.


## 3.2. Wavelet Interpolation Method for Squeeze Film Damping Equations

### 3.2.1. Construction of Basis Functions

In this subsection, the approximate solution ofPA(X,Y) is approximated by the following form:
(3.19)PA(X,Y)≈∑k=02J∑k′=02Jpk,k′JΘk,k′J(X,Y)=∑k=02J∑k′=02Jpk,k′Jθ̂J,k(X)θ̂J,k′(Y)=∑k=02J∑k′=02Jpk,k′Jθ̂(2JX-k)θ̂(2JY-k′),k,k′∈Z,
where the unknowns pk,k′J are the values of PA(X,Y) at the dyadic points X=k2-J, and  Y=k′2-J. The unknowns pk,k′J are complex.For the application of Galerkin method, (3.19) should be able to satisfy the boundary conditions. Substituting (3.19) into (3.18), leads to(3.20)∑k′=02Jp0,k′Jθ̂(0)θ̂(2JY-k′)=0⇒p0,k′J=0,fork′=0,1,2,…2J,∑k′=02Jp2J,k′Jθ̂(0)θ̂(2JY-k′)=0⇒p2J,k′J=0,fork′=0,1,2,…2J,∑k=12J-1pk,0Jθ̂(2JX-k)θ̂(0)=0⇒pk,0J=0,fork=1,2,…,(2J-1),∑k=12J-1pk,2JJθ̂(2JX-k)θ̂(0)=0⇒pk,2JJ=0,fork=1,2,…,(2J-1).
Thus (3.19) is rewritten as (3.21)PA(X,Y)=∑k=12J-1∑k′=12J-1pk,k′JΘk,k′J(X,Y)=∑k=12J-1∑k′=12J-1pk,k′Jθ̂(2JX-k)θ̂(2JY-k′).


### 3.2.2. Discretion of the Boundary Value Problem

The weak form functional of (3.17) is(3.22)W(PA)=∬Ω{12[(∂PA∂X)2+β2(∂PA∂Y)2+jσωPA2]-jσωPA}dXdY.
From the necessary conditions for the determination of the minimum W, we obtain(3.23)δW(PA)=∬Ω[∂δPA∂X∂PA∂X+β2∂δPA∂Y∂PA∂Y+jσωδPAPA-jσωδPA]dXdY=0.
Substituting (3.19) into (3.23), leads to(3.24)∑k=12J-1∑k′=12J-1{∬Ω[∂Θm,nJ∂X∂Θk,k′J∂X+β2∂Θm,nJ∂Y∂Θk,k′J∂Y+jσωΘm,nJΘk,k′J]dXdY}pk,k′J=jσω∬ΩΘm,nJdXdY,form,n=1,2,…,(2J-1).
This is a (2J-1)2×(2J-1)2 linear system(3.25)Θp=jσωE,
where p=[p1,1Jp1,2J⋯p1,2J-1Jp2,1Jp2,2J⋯p2J-1,2J-1J]T is an (2J-1)2×1 unknown coefficients’ vector, E=[∬ΩΘ1,1JdXdY∬ΩΘ1,2JdXdY⋯∬ΩΘ2J-1,2J-1JdXdY]T is a (2J-1)2×1 matrix, and Θ is a (2J-1)2×(2J-1)2 matrix. The entries in Θ are of the form(3.26)Θ̂(i,j)=Θ̂([(m-1)(2J-1)+n],[(k-1)(2J-1)+k′])=∬Ω[∂θ̂(2JX-m)∂Xθ̂(2JY-n)∂θ̂(2JX-k)∂Xθ̂(2JY-k′)+β2θ̂(2JX-m)∂θ̂(2JY-n)∂Yθ̂(2JX-k)∂θ̂(2JY-k′)∂Y+jσωθ̂(2JX-m)θ̂(2JY-n)θ̂(2JX-k)θ̂(2JY-k′)∂θ̂(2JX-m)∂X]dXdY.


### 3.2.3. Squeeze Film Damping of the Parallel Plate

The numerical solution of (3.25) can be written as(3.27)p=jσωΘ-1E.
The elements of p can be expressed as(3.28)pk,k′J=pk,k′J,R+jpk,k′J,Iform,n=1,2,…,(2J-1),
where pk,k′J,R and pk,k′J,I are the real and imaginary parts of pk,k′J, respectively. Using (3.21) and (3.28), the force acting on the plate owing to the pressure of the squeeze gas film can be rewritten as
(3.29)Pnon∬01(P̂-1)dXdY=AejωTS·Pnon∑k=12J-1∑k′=12J-1(pk,k′J,R+jpk,k′J,I)∬01θ̂(2JX-k)θ̂(2JY-k′)dXdY=Ka·Z(T)+Ca·dZ(T)dT,
where(3.30)Ka=Pnon∑k=12J-1∑k′=12J-1pk,k′J,R∬01θ̂(2JX-k)θ̂(2JY-k′)dXdY,Ca=PnonωS∑k=12J-1∑k′=12J-1pk,k′J,I∬01θ̂(2JX-k)θ̂(2JY-k′)dXdYKa·Z(T) and Ca·(dZ(T)/dT) are the spring and damping components of the force. Substituting (3.29), (3.14) and (3.15) into (3.8), we obtain(3.31)d2ZdT2+Ca·dZdT+(Ka+1)Z(T)=αV0v(T),Z(T)=AejωTS=aV0v0Ka+1·11-(ω2(S2/(Ka+1)))+(jω(CaS/(Ka+1)))ejωTS,
where S=mplate/kE=1/ωn. The quality factor and the damped natural frequency are expressed as(3.32)Qsqueeze=12ξ=Ka+1Ca,ωsqueeze=ωnKa+1.


## 3.2.1. Construction of Basis Functions

In this subsection, the approximate solution ofPA(X,Y) is approximated by the following form:
(3.19)PA(X,Y)≈∑k=02J∑k′=02Jpk,k′JΘk,k′J(X,Y)=∑k=02J∑k′=02Jpk,k′Jθ̂J,k(X)θ̂J,k′(Y)=∑k=02J∑k′=02Jpk,k′Jθ̂(2JX-k)θ̂(2JY-k′),k,k′∈Z,
where the unknowns pk,k′J are the values of PA(X,Y) at the dyadic points X=k2-J, and  Y=k′2-J. The unknowns pk,k′J are complex.For the application of Galerkin method, (3.19) should be able to satisfy the boundary conditions. Substituting (3.19) into (3.18), leads to(3.20)∑k′=02Jp0,k′Jθ̂(0)θ̂(2JY-k′)=0⇒p0,k′J=0,fork′=0,1,2,…2J,∑k′=02Jp2J,k′Jθ̂(0)θ̂(2JY-k′)=0⇒p2J,k′J=0,fork′=0,1,2,…2J,∑k=12J-1pk,0Jθ̂(2JX-k)θ̂(0)=0⇒pk,0J=0,fork=1,2,…,(2J-1),∑k=12J-1pk,2JJθ̂(2JX-k)θ̂(0)=0⇒pk,2JJ=0,fork=1,2,…,(2J-1).
Thus (3.19) is rewritten as (3.21)PA(X,Y)=∑k=12J-1∑k′=12J-1pk,k′JΘk,k′J(X,Y)=∑k=12J-1∑k′=12J-1pk,k′Jθ̂(2JX-k)θ̂(2JY-k′).


## 3.2.2. Discretion of the Boundary Value Problem

The weak form functional of (3.17) is(3.22)W(PA)=∬Ω{12[(∂PA∂X)2+β2(∂PA∂Y)2+jσωPA2]-jσωPA}dXdY.
From the necessary conditions for the determination of the minimum W, we obtain(3.23)δW(PA)=∬Ω[∂δPA∂X∂PA∂X+β2∂δPA∂Y∂PA∂Y+jσωδPAPA-jσωδPA]dXdY=0.
Substituting (3.19) into (3.23), leads to(3.24)∑k=12J-1∑k′=12J-1{∬Ω[∂Θm,nJ∂X∂Θk,k′J∂X+β2∂Θm,nJ∂Y∂Θk,k′J∂Y+jσωΘm,nJΘk,k′J]dXdY}pk,k′J=jσω∬ΩΘm,nJdXdY,form,n=1,2,…,(2J-1).
This is a (2J-1)2×(2J-1)2 linear system(3.25)Θp=jσωE,
where p=[p1,1Jp1,2J⋯p1,2J-1Jp2,1Jp2,2J⋯p2J-1,2J-1J]T is an (2J-1)2×1 unknown coefficients’ vector, E=[∬ΩΘ1,1JdXdY∬ΩΘ1,2JdXdY⋯∬ΩΘ2J-1,2J-1JdXdY]T is a (2J-1)2×1 matrix, and Θ is a (2J-1)2×(2J-1)2 matrix. The entries in Θ are of the form(3.26)Θ̂(i,j)=Θ̂([(m-1)(2J-1)+n],[(k-1)(2J-1)+k′])=∬Ω[∂θ̂(2JX-m)∂Xθ̂(2JY-n)∂θ̂(2JX-k)∂Xθ̂(2JY-k′)+β2θ̂(2JX-m)∂θ̂(2JY-n)∂Yθ̂(2JX-k)∂θ̂(2JY-k′)∂Y+jσωθ̂(2JX-m)θ̂(2JY-n)θ̂(2JX-k)θ̂(2JY-k′)∂θ̂(2JX-m)∂X]dXdY.


## 3.2.3. Squeeze Film Damping of the Parallel Plate

The numerical solution of (3.25) can be written as(3.27)p=jσωΘ-1E.
The elements of p can be expressed as(3.28)pk,k′J=pk,k′J,R+jpk,k′J,Iform,n=1,2,…,(2J-1),
where pk,k′J,R and pk,k′J,I are the real and imaginary parts of pk,k′J, respectively. Using (3.21) and (3.28), the force acting on the plate owing to the pressure of the squeeze gas film can be rewritten as
(3.29)Pnon∬01(P̂-1)dXdY=AejωTS·Pnon∑k=12J-1∑k′=12J-1(pk,k′J,R+jpk,k′J,I)∬01θ̂(2JX-k)θ̂(2JY-k′)dXdY=Ka·Z(T)+Ca·dZ(T)dT,
where(3.30)Ka=Pnon∑k=12J-1∑k′=12J-1pk,k′J,R∬01θ̂(2JX-k)θ̂(2JY-k′)dXdY,Ca=PnonωS∑k=12J-1∑k′=12J-1pk,k′J,I∬01θ̂(2JX-k)θ̂(2JY-k′)dXdYKa·Z(T) and Ca·(dZ(T)/dT) are the spring and damping components of the force. Substituting (3.29), (3.14) and (3.15) into (3.8), we obtain(3.31)d2ZdT2+Ca·dZdT+(Ka+1)Z(T)=αV0v(T),Z(T)=AejωTS=aV0v0Ka+1·11-(ω2(S2/(Ka+1)))+(jω(CaS/(Ka+1)))ejωTS,
where S=mplate/kE=1/ωn. The quality factor and the damped natural frequency are expressed as(3.32)Qsqueeze=12ξ=Ka+1Ca,ωsqueeze=ωnKa+1.


## 4. Wavelet Interpolation Galerkin Method for a Torsion Microplate    under the Effect of Squeeze Film Damping

A similar analysis as the one given for the parallel plate microresonator can be given for a torsion microplate.


### 4.1. Governing Equations

In this section, we examine the example of a rectangular torsion microplate under the effect of squeeze film damping. As shown in Figure2, the microplate is suspended by two torsion microbeams. lx, ly and hδ are the length, width and thickness of the plate. There are two pairs of electrodes between the microplate and the substrate. The locations of the two pairs of electrodes are symmetrical. x1 and x2 are the positions of the two pairs of electrodes. The thickness of the electrodes is neglected. On each pair of the electrodes, an equal dc voltage V0 and an equal ac voltage v(t) with opposite potential were applied. The rotation angle of the plate is composed of a static component to the dc voltage, denoted by γ0, and a small dynamic component due to the ac voltage, denoted by γ(t). In this case, γ0=0; thus the equation of the plate around V0 and γ0 can be written as
(4.1)Jγ̈+kT-Eγ=-(x22-x12)εlyV0g02v+∫-lx/2lx/2∫0l(p(x,y,t)-p0)xdydx,
where kT-E=[kT-2εlyV02((x23-x13)/3g03)] and kT is the stiffness of the two torsion microbeams. The pressure p(x,y,t) is governed by (3.5), where h(x,t)=g0+xγ(t). The pressure boundary conditions for the case in Figure 2 are(4.2)p(x,0,t)=p(x,ly,t)=p(-lx2,y,t)=p(lx2,y,t)=p0.A schematic drawing of a torsion microplate under the effect of squeeze film damping.
(a)
3-D view of the microplate(b)
Side view of the microplateFor convenience, we introduce the nondimensional variables(4.3)X=xlx+12,X1=x1lx+12,X2=x2lx+12,Y=yly,ϑ=γγmax,γmax=2g0lx,P̂=pp0,T=tS,H=hg0=1+2(X-12)ϑ,
where S=J/(kT-E)=1/ωn, ωn is the nature frequency of the plate. Substituting (4.3) into (4.1), (3.5) and (4.2), we obtain(4.4)ϑ̈+ϑ=-αV0v+Pnon∬01(P̂(X,Y,T)-1)(X-12)dXdY,(4.5)∂∂X(H3P̂∂P̂∂X)+β2∂∂Y(H3P̂∂P̂∂Y)=σS(H∂P̂∂T+P̂∂H∂T),
where α=(x22-x12)εly/kT-Eg02γmax,Pnon=p0lx2ly/kT-eγmax, σ=12ηlx2/g02p0, and  β=lx/ly. The nondimensional boundary conditions are(4.6)P̂(X,0,T)=P̂(X,1,T)=P̂(0,Y,T)=P̂(1,Y,T)=1.As mentioned above, the microplate is under small torsion oscillation aroundγ0=0 and therefore the pressure variation from ambient in the squeeze film is also small, P̂(X,Y,T) is given by
(4.7)P̂(X,Y,T)=pp0=1+P(X,Y,T),
where |P(X,Y,T)|≪1. Substituting (4.7) into (4.5), and linearizing the outcome around p0 and γ0, we obtain
(4.8)∂2P∂X2+β2∂2P∂Y2-σS∂P∂T=2σS(X-12)∂ϑ∂T.
The boundary conditions for the case are(4.9)P(X,0,T)=P(X,1,T)=P(0,Y,T)=P(1,Y,T)=0.For a harmonic excitation, the ac component voltagev(T) is given by(4.10)v(T)=v0ejωTS.
Correspondingly, the steady-state solution of (4.4) and (4.8) may be expressed by(4.11)ϑ(T)=AejωTS,P(X,Y,T)=A·PA(X,Y)ejωTS,
where A is the complex amplitude to be determined. Substituting (4.11) into (4.8), we obtain
(4.12)∂2PA(X,Y)∂X2+β2∂2PA(X,Y)∂Y2-jσωPA(X,Y)=j2σω(X-12).
The boundary conditions are(4.13)PA(X,0)=PA(X,1)=PA(0,Y)=PA(1,Y)=0.


### 4.2. Wavelet Interpolation Method for Squeeze Film Damping Equations

In this section, the approximate solution ofPA(X,Y) can be approximated by (3.21). The weak form functional of (4.12) is(4.14)W(PA)=∬Ω{12[(∂PA∂X)2+β2(∂PA∂Y)2+jσωPA2]+j2σω(X-12)PA}dXdY.
From the necessary conditions for the determination of the minimum W, we obtain
(4.15)δW(PA)=∬Ω[∂δPA∂X∂PA∂X+β2∂δPA∂Y∂PA∂Y+jσωδPAPA+j2σω(X-12)δPA]dXdY=0.
Substituting (3.21) into (4.15), leads to(4.16)∑k=12J-1∑k′=12J-1{∬Ω[∂Θm,nJ∂X∂Θk,k′J∂X+β2∂Θm,nJ∂Y∂Θk,k′J∂Y+jσωΘm,nJΘk,k′J]dXdY}pk,k′J=-j2σω∬Ω(X-12)Θm,nJdXdY,form,n=1,2,…,(2J-1).
This is a (2J-1)2×(2J-1)2 linear system(4.17)Θp=-jσωE,
where p=[p1,1Jp1,2J⋯p1,2J-1Jp2,1Jp2,2J⋯p2J-1,2J-1J]T is an (2J-1)2×1unknown coefficients’matrix, E=2[∬Ω(X-1/2)Θ1,1JdXdY∬Ω(X-1/2)Θ1,2JdXdY⋯∬Ω(X-(1/2))Θ2J-1,2J-1JdXdY]Tis a (2J-1)2×1 matrix, and Θ is a (2J-1)2×(2J-1)2 matrix. The entries in Θ are of the form(4.18)Θ̂(i,j)=Θ̂([(m-1)(2J-1)+n],[(k-1)(2J-1)+k′])=∬Ω[∂θ̂(2JX-m)∂Xθ̂(2JY-n)∂θ̂(2JX-k)∂Xθ̂(2JY-k′)+β2θ̂(2JX-m)∂θ̂(2JY-n)∂Yθ̂(2JX-k)∂θ̂(2JY-k′)∂Y+jσωθ̂(2JX-m)θ̂(2JY-n)θ̂(2JX-k)θ̂(2JY-k′)∂θ̂(2JX-m)∂X]dXdY.


#### 4.2.1. Squeeze Film Damping of the Torsion Plate

The numerical solution of (4.17) can be written as(4.19)p=-jσωΘ-1E.
The elements of p can be expressed as(4.20)pk,k′J=-pk,k′J,R-jpk,k′J,Iform,n=1,2,…,(2J-1),
where pk,k′J,R and pk,k′J,I are the real and imaginary parts of pk,k′J, respectively. Using (4.20) and (4.4), the force acting on the plate owing to the pressure of the squeeze gas film can be rewritten as(4.21)Pnon∬01(P̂-1)(X-12)dXdY=-Ka·ϑ(T)-Ca·dϑ(T)dT,
where Ka·ϑ(T) and Ca·(dϑ(T)/dT) are the spring and damping components of the force(4.22)Ka=Pnon∑k=12J-1∑k′=12J-1pk,k′J,R∬01θ̂(2JX-k)θ̂(2JY-k′)(X-12)dXdYCa=PnonωS∑k=12J-1∑k′=12J-1pk,k′J,I∬01θ̂(2JX-k)θ̂(2JY-k′)(X-12)dXdY
Substituting (4.21) into (4.4), leads to(4.23)d2ϑdT2+Ca·dϑdT+(Ka+1)ϑ=αV0v(T),ϑ(T)=AejωTS=aV0v0Ka+1·11-(ω2(S2/(Ka+1)))+(jω(CaS/(Ka+1)))ejωTS.
The quality factor and the damped natural frequency are given in(3.32)


## 4.1. Governing Equations

In this section, we examine the example of a rectangular torsion microplate under the effect of squeeze film damping. As shown in Figure2, the microplate is suspended by two torsion microbeams. lx, ly and hδ are the length, width and thickness of the plate. There are two pairs of electrodes between the microplate and the substrate. The locations of the two pairs of electrodes are symmetrical. x1 and x2 are the positions of the two pairs of electrodes. The thickness of the electrodes is neglected. On each pair of the electrodes, an equal dc voltage V0 and an equal ac voltage v(t) with opposite potential were applied. The rotation angle of the plate is composed of a static component to the dc voltage, denoted by γ0, and a small dynamic component due to the ac voltage, denoted by γ(t). In this case, γ0=0; thus the equation of the plate around V0 and γ0 can be written as
(4.1)Jγ̈+kT-Eγ=-(x22-x12)εlyV0g02v+∫-lx/2lx/2∫0l(p(x,y,t)-p0)xdydx,
where kT-E=[kT-2εlyV02((x23-x13)/3g03)] and kT is the stiffness of the two torsion microbeams. The pressure p(x,y,t) is governed by (3.5), where h(x,t)=g0+xγ(t). The pressure boundary conditions for the case in Figure 2 are(4.2)p(x,0,t)=p(x,ly,t)=p(-lx2,y,t)=p(lx2,y,t)=p0.A schematic drawing of a torsion microplate under the effect of squeeze film damping.
(a)
3-D view of the microplate(b)
Side view of the microplateFor convenience, we introduce the nondimensional variables(4.3)X=xlx+12,X1=x1lx+12,X2=x2lx+12,Y=yly,ϑ=γγmax,γmax=2g0lx,P̂=pp0,T=tS,H=hg0=1+2(X-12)ϑ,
where S=J/(kT-E)=1/ωn, ωn is the nature frequency of the plate. Substituting (4.3) into (4.1), (3.5) and (4.2), we obtain(4.4)ϑ̈+ϑ=-αV0v+Pnon∬01(P̂(X,Y,T)-1)(X-12)dXdY,(4.5)∂∂X(H3P̂∂P̂∂X)+β2∂∂Y(H3P̂∂P̂∂Y)=σS(H∂P̂∂T+P̂∂H∂T),
where α=(x22-x12)εly/kT-Eg02γmax,Pnon=p0lx2ly/kT-eγmax, σ=12ηlx2/g02p0, and  β=lx/ly. The nondimensional boundary conditions are(4.6)P̂(X,0,T)=P̂(X,1,T)=P̂(0,Y,T)=P̂(1,Y,T)=1.As mentioned above, the microplate is under small torsion oscillation aroundγ0=0 and therefore the pressure variation from ambient in the squeeze film is also small, P̂(X,Y,T) is given by
(4.7)P̂(X,Y,T)=pp0=1+P(X,Y,T),
where |P(X,Y,T)|≪1. Substituting (4.7) into (4.5), and linearizing the outcome around p0 and γ0, we obtain
(4.8)∂2P∂X2+β2∂2P∂Y2-σS∂P∂T=2σS(X-12)∂ϑ∂T.
The boundary conditions for the case are(4.9)P(X,0,T)=P(X,1,T)=P(0,Y,T)=P(1,Y,T)=0.For a harmonic excitation, the ac component voltagev(T) is given by(4.10)v(T)=v0ejωTS.
Correspondingly, the steady-state solution of (4.4) and (4.8) may be expressed by(4.11)ϑ(T)=AejωTS,P(X,Y,T)=A·PA(X,Y)ejωTS,
where A is the complex amplitude to be determined. Substituting (4.11) into (4.8), we obtain
(4.12)∂2PA(X,Y)∂X2+β2∂2PA(X,Y)∂Y2-jσωPA(X,Y)=j2σω(X-12).
The boundary conditions are(4.13)PA(X,0)=PA(X,1)=PA(0,Y)=PA(1,Y)=0.


## 4.2. Wavelet Interpolation Method for Squeeze Film Damping Equations

In this section, the approximate solution ofPA(X,Y) can be approximated by (3.21). The weak form functional of (4.12) is(4.14)W(PA)=∬Ω{12[(∂PA∂X)2+β2(∂PA∂Y)2+jσωPA2]+j2σω(X-12)PA}dXdY.
From the necessary conditions for the determination of the minimum W, we obtain
(4.15)δW(PA)=∬Ω[∂δPA∂X∂PA∂X+β2∂δPA∂Y∂PA∂Y+jσωδPAPA+j2σω(X-12)δPA]dXdY=0.
Substituting (3.21) into (4.15), leads to(4.16)∑k=12J-1∑k′=12J-1{∬Ω[∂Θm,nJ∂X∂Θk,k′J∂X+β2∂Θm,nJ∂Y∂Θk,k′J∂Y+jσωΘm,nJΘk,k′J]dXdY}pk,k′J=-j2σω∬Ω(X-12)Θm,nJdXdY,form,n=1,2,…,(2J-1).
This is a (2J-1)2×(2J-1)2 linear system(4.17)Θp=-jσωE,
where p=[p1,1Jp1,2J⋯p1,2J-1Jp2,1Jp2,2J⋯p2J-1,2J-1J]T is an (2J-1)2×1unknown coefficients’matrix, E=2[∬Ω(X-1/2)Θ1,1JdXdY∬Ω(X-1/2)Θ1,2JdXdY⋯∬Ω(X-(1/2))Θ2J-1,2J-1JdXdY]Tis a (2J-1)2×1 matrix, and Θ is a (2J-1)2×(2J-1)2 matrix. The entries in Θ are of the form(4.18)Θ̂(i,j)=Θ̂([(m-1)(2J-1)+n],[(k-1)(2J-1)+k′])=∬Ω[∂θ̂(2JX-m)∂Xθ̂(2JY-n)∂θ̂(2JX-k)∂Xθ̂(2JY-k′)+β2θ̂(2JX-m)∂θ̂(2JY-n)∂Yθ̂(2JX-k)∂θ̂(2JY-k′)∂Y+jσωθ̂(2JX-m)θ̂(2JY-n)θ̂(2JX-k)θ̂(2JY-k′)∂θ̂(2JX-m)∂X]dXdY.


### 4.2.1. Squeeze Film Damping of the Torsion Plate

The numerical solution of (4.17) can be written as(4.19)p=-jσωΘ-1E.
The elements of p can be expressed as(4.20)pk,k′J=-pk,k′J,R-jpk,k′J,Iform,n=1,2,…,(2J-1),
where pk,k′J,R and pk,k′J,I are the real and imaginary parts of pk,k′J, respectively. Using (4.20) and (4.4), the force acting on the plate owing to the pressure of the squeeze gas film can be rewritten as(4.21)Pnon∬01(P̂-1)(X-12)dXdY=-Ka·ϑ(T)-Ca·dϑ(T)dT,
where Ka·ϑ(T) and Ca·(dϑ(T)/dT) are the spring and damping components of the force(4.22)Ka=Pnon∑k=12J-1∑k′=12J-1pk,k′J,R∬01θ̂(2JX-k)θ̂(2JY-k′)(X-12)dXdYCa=PnonωS∑k=12J-1∑k′=12J-1pk,k′J,I∬01θ̂(2JX-k)θ̂(2JY-k′)(X-12)dXdY
Substituting (4.21) into (4.4), leads to(4.23)d2ϑdT2+Ca·dϑdT+(Ka+1)ϑ=αV0v(T),ϑ(T)=AejωTS=aV0v0Ka+1·11-(ω2(S2/(Ka+1)))+(jω(CaS/(Ka+1)))ejωTS.
The quality factor and the damped natural frequency are given in(3.32)


## 4.2.1. Squeeze Film Damping of the Torsion Plate

The numerical solution of (4.17) can be written as(4.19)p=-jσωΘ-1E.
The elements of p can be expressed as(4.20)pk,k′J=-pk,k′J,R-jpk,k′J,Iform,n=1,2,…,(2J-1),
where pk,k′J,R and pk,k′J,I are the real and imaginary parts of pk,k′J, respectively. Using (4.20) and (4.4), the force acting on the plate owing to the pressure of the squeeze gas film can be rewritten as(4.21)Pnon∬01(P̂-1)(X-12)dXdY=-Ka·ϑ(T)-Ca·dϑ(T)dT,
where Ka·ϑ(T) and Ca·(dϑ(T)/dT) are the spring and damping components of the force(4.22)Ka=Pnon∑k=12J-1∑k′=12J-1pk,k′J,R∬01θ̂(2JX-k)θ̂(2JY-k′)(X-12)dXdYCa=PnonωS∑k=12J-1∑k′=12J-1pk,k′J,I∬01θ̂(2JX-k)θ̂(2JY-k′)(X-12)dXdY
Substituting (4.21) into (4.4), leads to(4.23)d2ϑdT2+Ca·dϑdT+(Ka+1)ϑ=αV0v(T),ϑ(T)=AejωTS=aV0v0Ka+1·11-(ω2(S2/(Ka+1)))+(jω(CaS/(Ka+1)))ejωTS.
The quality factor and the damped natural frequency are given in(3.32)


## 5. Comparisons with Experiments

Veijola et al. [26] conducted experiments to measure the frequency response of an accelerometer under the effect of squeeze film damping. Minikes et al. [27] measured the quality factors of two torsion rectangular mirrors at low pressure. In this section, the experimental results presented by Veijola et al. [26] and Minikes et al. [27] were used to verify the wavelet interpolation Galerkin method.


### 5.1. Comparisons with the Experimental Results of Veijola et al. [26]

In  [26], Veijola et al. simulated the frequency response of an accelerometer with a spring-mass- damper model with a parallel-plate electrostatic force. The damping coefficient was estimated by the Blech model [28]. The spring constants and the gas pressures were estimated by curve fitting the experimental measurements. They compared their simulations with experimental data and found good agreement. The parameters for the accelerometer are listed in Table 1.Table 1
The parameters of the accelerometer presented by Veijola et al. [26].
ParametersValuesMassmplate4.9×10-6 kgSpring constantkspring212.1 N/mGap spacingg03.95μmLength of the moving massslx2 960μmWidth of the moving masssly1 780μmAmbient pressurep011 PaBias voltageV09 VEffective viscosityηeff10.2×10-9 N·s·m-2In this subsection, we use the wavelet interpolation Galerkin method to predict the frequency response of the accelerometer. Various numerical tests have been conducted by changing the degree of the Daubechies waveletL and the number J of the scale. Better accuracy can be achieved by increasing L and J. The higher L is, the smoother the scaling function becomes. The price for the high smoothness is that its supporting domain gets larger. The higher J is, the more accurate the solution becomes. The number of differential equations and the CPU time increase significantly as J increases. In this work, only the solutions for L=6 and J=4 are presented, as the results for higher resolutions are indistinguishable from the exact solution.For comparison purpose, we give the frequency responses of the accelerometer calculated by Blech’s model [28] and the finite difference method, respectively. Blech [28] expanded the air film pressure into an assumed double sine series and derived an analytical expression for the spring and damping forces. In this subsection, the number of terms for the double sine series is taken as (5,5), which shows good convergence. For the finite difference method, we use the following approximate formulae for a node (i,j) on the microplate:
(5.1)∂PA∂X|i,j=PA(i+1,j)-PA(i-1,j)2ΔX,∂PA∂Y|i,j=PA(i,j+1)-PA(i,j-1)2ΔY,∂2PA∂X2|i,j=PA(i+1,j)-2PA(i,j)+PA(i-1,j)ΔX2,∂2PA∂Y2|i,j=PA(i,j+1)-2PA(i,j)+PA(i,j-1)ΔY2.
In this subsection, we assume that ΔX=ΔY=1/2J=1/26; thus the element size of the finite difference method is equal to the wavelet interpolation Galerkin method. Substituting (5.1) into (3.17), we obtain(5.2)PA(i+1,j)+PA(i-1,j)+β2PA(i,j+1)+β2PA(i,j-1)ΔX2-(2+2β2ΔX2+jσω)PA(i,j)=-jσω.Figure3 shows the comparisons of the frequency response obtained by different methods. As expected, the wavelet interpolation Galerkin method, Blech’s model and the finite difference method give almost same results. The three results agree well with the experimental results [26] except for one data at resonance peak of the amplitude frequency response. The reason for this discrepancy is that the damping coefficient is slightly underestimated by the three methods, respectively. Table 2 shows the Comparison of the damping obtained by different methods. In the experiment [26], the squeeze film damping is dominant. Obviously, the accuracy of the finite difference method is less than the wavelet interpolation Galerkin method and Blech’s model. The wavelet interpolation Galerkin method and the Blech model give almost identical results. Figure 4 shows the real part and the imaginary part of the air film pressure calculated by the wavelet interpolation Galerkin method.Table 2
A Comparison of the damping ratios obtained by different methods with the experimental data [26].
MethodsDamping ratioPeak value (dB)Experimental data [26]0.023918.3Blech’s model (error)0.0159 (33.5%)21.7The wavelet interpolation Galerkin method (error)0.0155 (35.1%)21.9The finite difference method (error)0.0136 (43.1%)23.2Comparisons of the frequency response obtained by the wavelet Galerkin method, the Blech model and the experimental data of Veijola et al. [26].
(a)
Amplitude response(b)
Phase responseThe air film pressure distribution calculated by the wavelet interpolation Galerkin method.
(a)
The real part of the pressure distribution(b)
The imaginary part of the pressure distribution


### 5.2. Comparisons with the Experimental Results of Minikes et al. [27]

Minikes et al. [27] measured the quality factors of two rectangular torsion mirrors at low pressure and plotted the curves of the quality factors as a function of air pressure in the range from 10−2 torr to 102 torr. The structure of the two torsion mirrors is identical with the structure shown in Figure 2. The two torsion mirrors have similar dimensions in terms of surface area and inertial moment, but have different gaps between the mirror and the actuation electrodes. The parameters for the two torsion mirrors are listed in Table 3.Table 3
The parameters of two torsion mirrors [27].
Mirror 1Mirror 2Width of the mirrorly500μm500μmLength of the mirrorlx500μm500μmThickness of the mirrorhδ30μm30μmDensity of the mirrorρ2300 kg/m32300 kg/m3Gap spacingg028μm13μmTorsional natural frequencyfn13092.56 Hz12824.87 HzIn Table3, the values of the two torsional natural frequencies are measured under the dc bias voltage. Based on the two torsional natural frequencies and two moments of inertia, we determined  the torsional stiffness kT-E. The extracted torsional stiffness kT-E for mirror 1 and mirror 2 are 2.461×10−6 and 2.362×10−6 N·m/rad, respectively. Now we use the wavelet interpolation Galerkin method to predict the quality factors of the two torsion mirrors. Various numerical tests have been conducted by changing the degree of the Daubechies wavelet L and the number J of the scale. Only the solutions for L=6 and J=4 are presented, as the results for higher resolutions are indistinguishable from the exact solution.For comparison purpose, we give the quality factors calculated by Pan’s model [15] and the finite difference method, respectively. Pan et al. [15] expanded the air film pressure into an assumed double sine series and derived an analytical expression for the spring and damping torques. In this work, the number of terms for the double sine series is taken as (6,5), which shows good convergence. For the finite difference method, we use (5.1) and (4.12) for a node (i,j) on the mirror; thus obtain(5.3)PA(i+1,j)+PA(i-1,j)+β2PA(i,j+1)+β2PA(i,j-1)ΔX2-(2+2β2ΔX2+jσω)PA(i,j)=j2σω(Xi-12).
In this subsection, we assume that ΔX=ΔY=1/2J =1/26; thus the finite difference method yields the same grids as the wavelet interpolation Galerkin method.Tables4 and 5 show the comparisons of quality factors obtained by the wavelet interpolation Galerkin method, the finite difference method and Pan’s model for the mirror 1 and 2, respectively. As shown in Tables 4 and 5, the qualify factors obtained by the wavelet interpolation Galerkin method are almost 1%~4% higher than Pan’s model. However the qualify factors obtained by the finite difference method are almost 10%~15% higher than Pan’s model. The result of the wavelet interpolation Galerkin method matches the result of Pan’s model better than that of the finite difference method.Table 4
A comparison of quality factors obtained by Pan’s model, the wavelet interpolation Galerkin method and the finite difference method for mirror 1.
p0(torr)Pan’s modelThe wavelet interpolation Galerkin method (error)The finite difference method (error)0.085.784×1045.987× 104 (3.5%)6.586 × 104 (13.9%)0.104.483×1044.627× 104 (3.2%)5.160 × 104 (15.1%)0.507.070×1037.330× 103 (3.7%)8.006 × 103 (13.2%)13.306×1033.388 × 103 (2.5%)3.694 × 103 (11.7%)31.061×1031.080 × 103 (1.8%)1.173 × 103 (10.6%)65.680×1025.895 × 102 (3.8%)6.526 × 102 (14.9%)103.979×1024.052 × 102 (1.8%)4.454 × 102 (11.9%)302.451×1022.505 × 102 (2.2%)2.721 × 102 (11.0%)602.074 × 1022.160 × 102 (4.1%)2.327 × 102 (12.2%)1001.994 × 1022.035 × 102 (2.1%)2.240 × 102 (12.3%)7601.850 × 1021.896 × 102 (2.5%)2.109 × 102  (14.0%)Table 5
A comparison of quality factors obtained by Pan’s model, the wavelet interpolation Galerkin method and the finite difference method for mirror 2.
p0(torr)Pan’s modelThe wavelet interpolation Galerkin method (error)The finite difference method (error)0.081.645×1041.679× 104 (2.1%)1.847 × 104 (12.3%)0.101.268 × 1041.297 × 104 (2.3%)1.458 × 104 (15.0%)0.502.135 × 1032.208 × 103 (3.4%)2.352 × 103 (10.2%)18.919 × 1029.195 × 102 (3.1%)1.010 × 103 (13.2%)32.681 × 1022.712 × 102 (1.2%)3.012 × 102 (12.3%)61.282 × 1021.318 × 102 (2.8%)1.463 × 102 (14.1%)108.004 × 1018.118 × 101 (1.4%)9.092 × 101 (13.6%)303.461 × 1013.599 × 101 (4.0%)3.959 × 101 (14.4%)602.574 × 1012.628 × 101 (2.1%)2.970 × 101 (15.4%)1002.255 × 1012.278 × 101 (1.0%)2.593 × 101 (15.0%)7601.820 × 1011.900 × 101 (4.4%)2.055 × 101 (12.9%)Figure5 shows the comparisons of quality factors obtained by different methods. As expected, the wavelet interpolation Galerkin method, the finite difference method and Pan’s model give almost same results. Above p0 = 10 torr, the viscous damping is dominant, the three methods give results in good agreement with the experimental results [27]. Below p0 = 10 torr, the accuracy of the three methods decreases as the pressure decreases. The main reason for this trend are as follows. The three methods are based on Reynolds equation, which is derived from the Navier-Stokes equations and the continuity equation. The main assumption is that the gas in the gap can be treated as a continuum. Below p0 = 10 torr, the Knudsen number Kn>0.1, the gas in the gap cannot be treated as a continuum. Thus the three methods fail to give a good prediction. Below p0 = 0.1 torr, the influence of squeeze film damping begins to vanish and the qualify factors reaches a plateau that is dominated by the intrinsic damping.Comparison of the quality factors obtained by the wavelet Galerkin method, the Pan model and the experimental data of Minikes et al. [27].
(a)
Mirror 1 (g0= 28 μm)(b)
Mirror 2 (g0= 13 μm)Tables6 and 7 list the comparisons of quality factors between p = 10 and 760 torr for mirrors 1 and 2, respectively. Obviously, the accuracy of the finite difference method is less than the wavelet interpolation Galerkin method and Pan’s model. Figure 6 shows the real part and the imaginary part of the air film pressure calculated by the wavelet interpolation Galerkin method at p0 = 1 torr. The air film pressure looks similar to the results calculated by Pan’s model.Table 6
A comparison of quality factors betweenp0 = 10 and 760 torr for mirror 1.
p0 (torr)Experiment data [27]Pan’s model (error)The wavelet interpolation Galerkin method (error)The finite difference method (error)10345398 (15.4%)405 (17.4%)445 (29.0%)20275289 (5.1%)301 (9.4%)311 (13.1%)40229236 (3.1%)241 (5.2%)252 (10.0%)760179185 (3.4%)190 (6.1%)211 (17.9%)Table 7
A comparison of quality factors betweenp0 = 10 and 760 torr for mirror 2.
p0(torr)Expertiment data [27]Pan’s model (error)The wavelet interpolation Galerkin method (error)The finite difference method (error)106680.0 (21.2%)81.2 (23.0%)90.9 (37.7%)204048.8 (22.0%)50.1 (25.3%)54.0 (35.0%)402730.1 (11.5%)31.0 (14.8%)33.5 (24.1%)7601818.2 (1.1%)19.0 (5.6%)20.6 (14.4%)The air film pressure distribution calculated by the wavelet interpolation Galerkin method atp0 = 1 torr.
(a)
The real part of mirror 1(b)
The imaginary part of mirror 1(c)
The real part of mirror 2(d)
The imaginary part of mirror 2


## 5.1. Comparisons with the Experimental Results of Veijola et al. [26]

In  [26], Veijola et al. simulated the frequency response of an accelerometer with a spring-mass- damper model with a parallel-plate electrostatic force. The damping coefficient was estimated by the Blech model [28]. The spring constants and the gas pressures were estimated by curve fitting the experimental measurements. They compared their simulations with experimental data and found good agreement. The parameters for the accelerometer are listed in Table 1.Table 1
The parameters of the accelerometer presented by Veijola et al. [26].
ParametersValuesMassmplate4.9×10-6 kgSpring constantkspring212.1 N/mGap spacingg03.95μmLength of the moving massslx2 960μmWidth of the moving masssly1 780μmAmbient pressurep011 PaBias voltageV09 VEffective viscosityηeff10.2×10-9 N·s·m-2In this subsection, we use the wavelet interpolation Galerkin method to predict the frequency response of the accelerometer. Various numerical tests have been conducted by changing the degree of the Daubechies waveletL and the number J of the scale. Better accuracy can be achieved by increasing L and J. The higher L is, the smoother the scaling function becomes. The price for the high smoothness is that its supporting domain gets larger. The higher J is, the more accurate the solution becomes. The number of differential equations and the CPU time increase significantly as J increases. In this work, only the solutions for L=6 and J=4 are presented, as the results for higher resolutions are indistinguishable from the exact solution.For comparison purpose, we give the frequency responses of the accelerometer calculated by Blech’s model [28] and the finite difference method, respectively. Blech [28] expanded the air film pressure into an assumed double sine series and derived an analytical expression for the spring and damping forces. In this subsection, the number of terms for the double sine series is taken as (5,5), which shows good convergence. For the finite difference method, we use the following approximate formulae for a node (i,j) on the microplate:
(5.1)∂PA∂X|i,j=PA(i+1,j)-PA(i-1,j)2ΔX,∂PA∂Y|i,j=PA(i,j+1)-PA(i,j-1)2ΔY,∂2PA∂X2|i,j=PA(i+1,j)-2PA(i,j)+PA(i-1,j)ΔX2,∂2PA∂Y2|i,j=PA(i,j+1)-2PA(i,j)+PA(i,j-1)ΔY2.
In this subsection, we assume that ΔX=ΔY=1/2J=1/26; thus the element size of the finite difference method is equal to the wavelet interpolation Galerkin method. Substituting (5.1) into (3.17), we obtain(5.2)PA(i+1,j)+PA(i-1,j)+β2PA(i,j+1)+β2PA(i,j-1)ΔX2-(2+2β2ΔX2+jσω)PA(i,j)=-jσω.Figure3 shows the comparisons of the frequency response obtained by different methods. As expected, the wavelet interpolation Galerkin method, Blech’s model and the finite difference method give almost same results. The three results agree well with the experimental results [26] except for one data at resonance peak of the amplitude frequency response. The reason for this discrepancy is that the damping coefficient is slightly underestimated by the three methods, respectively. Table 2 shows the Comparison of the damping obtained by different methods. In the experiment [26], the squeeze film damping is dominant. Obviously, the accuracy of the finite difference method is less than the wavelet interpolation Galerkin method and Blech’s model. The wavelet interpolation Galerkin method and the Blech model give almost identical results. Figure 4 shows the real part and the imaginary part of the air film pressure calculated by the wavelet interpolation Galerkin method.Table 2
A Comparison of the damping ratios obtained by different methods with the experimental data [26].
MethodsDamping ratioPeak value (dB)Experimental data [26]0.023918.3Blech’s model (error)0.0159 (33.5%)21.7The wavelet interpolation Galerkin method (error)0.0155 (35.1%)21.9The finite difference method (error)0.0136 (43.1%)23.2Comparisons of the frequency response obtained by the wavelet Galerkin method, the Blech model and the experimental data of Veijola et al. [26].
(a)
Amplitude response(b)
Phase responseThe air film pressure distribution calculated by the wavelet interpolation Galerkin method.
(a)
The real part of the pressure distribution(b)
The imaginary part of the pressure distribution


## 5.2. Comparisons with the Experimental Results of Minikes et al. [27]

Minikes et al. [27] measured the quality factors of two rectangular torsion mirrors at low pressure and plotted the curves of the quality factors as a function of air pressure in the range from 10−2 torr to 102 torr. The structure of the two torsion mirrors is identical with the structure shown in Figure 2. The two torsion mirrors have similar dimensions in terms of surface area and inertial moment, but have different gaps between the mirror and the actuation electrodes. The parameters for the two torsion mirrors are listed in Table 3.Table 3
The parameters of two torsion mirrors [27].
Mirror 1Mirror 2Width of the mirrorly500μm500μmLength of the mirrorlx500μm500μmThickness of the mirrorhδ30μm30μmDensity of the mirrorρ2300 kg/m32300 kg/m3Gap spacingg028μm13μmTorsional natural frequencyfn13092.56 Hz12824.87 HzIn Table3, the values of the two torsional natural frequencies are measured under the dc bias voltage. Based on the two torsional natural frequencies and two moments of inertia, we determined  the torsional stiffness kT-E. The extracted torsional stiffness kT-E for mirror 1 and mirror 2 are 2.461×10−6 and 2.362×10−6 N·m/rad, respectively. Now we use the wavelet interpolation Galerkin method to predict the quality factors of the two torsion mirrors. Various numerical tests have been conducted by changing the degree of the Daubechies wavelet L and the number J of the scale. Only the solutions for L=6 and J=4 are presented, as the results for higher resolutions are indistinguishable from the exact solution.For comparison purpose, we give the quality factors calculated by Pan’s model [15] and the finite difference method, respectively. Pan et al. [15] expanded the air film pressure into an assumed double sine series and derived an analytical expression for the spring and damping torques. In this work, the number of terms for the double sine series is taken as (6,5), which shows good convergence. For the finite difference method, we use (5.1) and (4.12) for a node (i,j) on the mirror; thus obtain(5.3)PA(i+1,j)+PA(i-1,j)+β2PA(i,j+1)+β2PA(i,j-1)ΔX2-(2+2β2ΔX2+jσω)PA(i,j)=j2σω(Xi-12).
In this subsection, we assume that ΔX=ΔY=1/2J =1/26; thus the finite difference method yields the same grids as the wavelet interpolation Galerkin method.Tables4 and 5 show the comparisons of quality factors obtained by the wavelet interpolation Galerkin method, the finite difference method and Pan’s model for the mirror 1 and 2, respectively. As shown in Tables 4 and 5, the qualify factors obtained by the wavelet interpolation Galerkin method are almost 1%~4% higher than Pan’s model. However the qualify factors obtained by the finite difference method are almost 10%~15% higher than Pan’s model. The result of the wavelet interpolation Galerkin method matches the result of Pan’s model better than that of the finite difference method.Table 4
A comparison of quality factors obtained by Pan’s model, the wavelet interpolation Galerkin method and the finite difference method for mirror 1.
p0(torr)Pan’s modelThe wavelet interpolation Galerkin method (error)The finite difference method (error)0.085.784×1045.987× 104 (3.5%)6.586 × 104 (13.9%)0.104.483×1044.627× 104 (3.2%)5.160 × 104 (15.1%)0.507.070×1037.330× 103 (3.7%)8.006 × 103 (13.2%)13.306×1033.388 × 103 (2.5%)3.694 × 103 (11.7%)31.061×1031.080 × 103 (1.8%)1.173 × 103 (10.6%)65.680×1025.895 × 102 (3.8%)6.526 × 102 (14.9%)103.979×1024.052 × 102 (1.8%)4.454 × 102 (11.9%)302.451×1022.505 × 102 (2.2%)2.721 × 102 (11.0%)602.074 × 1022.160 × 102 (4.1%)2.327 × 102 (12.2%)1001.994 × 1022.035 × 102 (2.1%)2.240 × 102 (12.3%)7601.850 × 1021.896 × 102 (2.5%)2.109 × 102  (14.0%)Table 5
A comparison of quality factors obtained by Pan’s model, the wavelet interpolation Galerkin method and the finite difference method for mirror 2.
p0(torr)Pan’s modelThe wavelet interpolation Galerkin method (error)The finite difference method (error)0.081.645×1041.679× 104 (2.1%)1.847 × 104 (12.3%)0.101.268 × 1041.297 × 104 (2.3%)1.458 × 104 (15.0%)0.502.135 × 1032.208 × 103 (3.4%)2.352 × 103 (10.2%)18.919 × 1029.195 × 102 (3.1%)1.010 × 103 (13.2%)32.681 × 1022.712 × 102 (1.2%)3.012 × 102 (12.3%)61.282 × 1021.318 × 102 (2.8%)1.463 × 102 (14.1%)108.004 × 1018.118 × 101 (1.4%)9.092 × 101 (13.6%)303.461 × 1013.599 × 101 (4.0%)3.959 × 101 (14.4%)602.574 × 1012.628 × 101 (2.1%)2.970 × 101 (15.4%)1002.255 × 1012.278 × 101 (1.0%)2.593 × 101 (15.0%)7601.820 × 1011.900 × 101 (4.4%)2.055 × 101 (12.9%)Figure5 shows the comparisons of quality factors obtained by different methods. As expected, the wavelet interpolation Galerkin method, the finite difference method and Pan’s model give almost same results. Above p0 = 10 torr, the viscous damping is dominant, the three methods give results in good agreement with the experimental results [27]. Below p0 = 10 torr, the accuracy of the three methods decreases as the pressure decreases. The main reason for this trend are as follows. The three methods are based on Reynolds equation, which is derived from the Navier-Stokes equations and the continuity equation. The main assumption is that the gas in the gap can be treated as a continuum. Below p0 = 10 torr, the Knudsen number Kn>0.1, the gas in the gap cannot be treated as a continuum. Thus the three methods fail to give a good prediction. Below p0 = 0.1 torr, the influence of squeeze film damping begins to vanish and the qualify factors reaches a plateau that is dominated by the intrinsic damping.Comparison of the quality factors obtained by the wavelet Galerkin method, the Pan model and the experimental data of Minikes et al. [27].
(a)
Mirror 1 (g0= 28 μm)(b)
Mirror 2 (g0= 13 μm)Tables6 and 7 list the comparisons of quality factors between p = 10 and 760 torr for mirrors 1 and 2, respectively. Obviously, the accuracy of the finite difference method is less than the wavelet interpolation Galerkin method and Pan’s model. Figure 6 shows the real part and the imaginary part of the air film pressure calculated by the wavelet interpolation Galerkin method at p0 = 1 torr. The air film pressure looks similar to the results calculated by Pan’s model.Table 6
A comparison of quality factors betweenp0 = 10 and 760 torr for mirror 1.
p0 (torr)Experiment data [27]Pan’s model (error)The wavelet interpolation Galerkin method (error)The finite difference method (error)10345398 (15.4%)405 (17.4%)445 (29.0%)20275289 (5.1%)301 (9.4%)311 (13.1%)40229236 (3.1%)241 (5.2%)252 (10.0%)760179185 (3.4%)190 (6.1%)211 (17.9%)Table 7
A comparison of quality factors betweenp0 = 10 and 760 torr for mirror 2.
p0(torr)Expertiment data [27]Pan’s model (error)The wavelet interpolation Galerkin method (error)The finite difference method (error)106680.0 (21.2%)81.2 (23.0%)90.9 (37.7%)204048.8 (22.0%)50.1 (25.3%)54.0 (35.0%)402730.1 (11.5%)31.0 (14.8%)33.5 (24.1%)7601818.2 (1.1%)19.0 (5.6%)20.6 (14.4%)The air film pressure distribution calculated by the wavelet interpolation Galerkin method atp0 = 1 torr.
(a)
The real part of mirror 1(b)
The imaginary part of mirror 1(c)
The real part of mirror 2(d)
The imaginary part of mirror 2


## 6. Summary and Conclusions

A new wavelet interpolation Galerkin method has been developed for the numerical simulation of MEMS devices under the effect of squeeze film damping. The air film pressure are expressed as a linear combination of a class of interpolating functions generated by autocorrelation of the usual compactly supported Daubechies scaling functions. To the best of our knowledge, this is the first time that wavelets have been used as basis functions for solving the PDEs of MEMS devices. As opposed to the previous wavelet-based methods that are all limited in one energy domain, the MEMS devices in the paper involve two coupled energy domains. Two typical electrically actuated micro devices with squeeze film damping effect are examined respectively to illustrate the wavelet interpolation Galerkin method. The method is validated by comparing its results with available theoretical and experimental results. The accuracy of the method is higher than the finite difference method.In this paper, the wavelet interpolation Galerkin method is not suitable to solve problems defined on nonrectangular domains, since higher-dimensional wavelets are constructed by employing the tensor product of the one-dimensional wavelets and so their application is restricted to rectangular domains. In this paper, both trial and weight functions are a class of interpolating functions generated by autocorrelation of the first-generation wavelets. Future area of research is based on the second-generation wavelets [29]. The main advantage of the second-generation wavelets is that the wavelets are constructed in the spatial domain and can be custom designed for complex domains. This work is currently under way.



---

*Source: 586718-2010-01-21.xml*