- Email: [email protected]
A variation of local point interpolation method (vLPIM) for analysis of microelectromechanical systems (MEMS) device
Engineering Analysis with Boundary Elements 28 (2004) 1261–1270 www.elsevier.com/locate/enganabound
A variation of local point interpolation method (v LPIM) for analysis of microelectromechanical systems (MEMS) device Hua Lia,*, Q.X. Wanga, K.Y. Lama,b a
Institute of High Performance Computing, National University of Singapore, 1 Science Park Road, No. 01-01 The Capricorn, Singapore Science Park II, Singapore, Singapore 117528 b Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore, Singapore 119260 Received 28 May 2003; revised 18 August 2003; accepted 18 August 2003 Available online 16 March 2004
Abstract In this paper, a variation of local point interpolation method (v LPIM) is presented for meshless analysis of microelectromechanical systems (MEMS) devices, in which the weight function is assigned by high-order spline functions, instead of the Galerkin formulation of the local point interpolation method (LPIM). The present v LPIM with Timoshenko beam theory are then used to analyze MEMS. In the conventional MEMS design and simulation, time-consuming mesh generation is required. Recently, a number of meshless techniques have been applied to MEMS simulations since they do not require mesh generation. As one of them, the LPIM is based on the local weak form of the partial differential equations and easily implements the boundary conditions due to the shape functions constructed by the PIM with the delta function properties. Several MEMS devices are studied by the present v LPIM. The simulated results are compared with those by other simulation approaches, and it is shown that the present meshless method is very efficient and accurate for the analysis of MEMS devices. q 2004 Elsevier Ltd. All rights reserved. Keywords: Microelectromechanical system; Meshless method; Point interpolation; Static and dynamic analyses
1. Introduction In the past years, attention to microelectromechanical systems (MEMS) has a phenomenal growth because of its wide application in automotive systems (transducers and accelerometers), avionics and aerospace (microscale actuators and sensors), manufacturing and fabrication (microsmart robots), medicine and bioengineering (DNA and genetic code analysis and synthesis, drug delivery, diagnostics and imaging), etc. [1]. However, the development of the increasingly complicated MEMS demands sophisticated simulation techniques [2,3]. The traditional finite element methods or computer-aided design (CAD) systems require the generation of the threedimensional mesh, since the MEMS devices typically involve multiple coupled energy domains and media. It is obvious that mesh generation is computationally expensive and mesh refinement is very difficult, especially for the MEMS devices with complicated geometries and * Corresponding author. Tel.: þ 65-64191-561; fax: þ 65-64191-380. E-mail address: [email protected] (H. Li). 0955-7997/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2003.08.008
multiphysic domains. Therefore, it is necessary to find alternatives to the traditional methods. With the development and the application of the meshless approaches, a lot of meshless modelings and simulations for MEMS (such as microswitches, microsensors/microactuators, etc.) [4,5] have been carried out. Generally, the meshless methods are grouped into two categories according to using or not using integration [6]. The first category is the meshless methods, which do not require integration, based on the strong form of the partial differential equations (PDEs). It includes the smooth particle hydrodynamics [7], the generalized finite difference method [8], etc. The second category is the meshless methods based on the weak form of PDEs, including the element-free Garlerkin method [9], the point interpolation method (PIM) [10], and other methods referred to Refs. [6,11]. Furthermore, in order to avoid the global integration, the local weak forms are also used for the development of local meshless techniques, such as the meshless Petrov – Galerkin method (MLPG) [12 – 15] and the family of local point interpolation method (LPIM) [6,16,17].
1262
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
In the application of the local meshless methods, the choice of the trial function and test function is very important to ensure the computational simplicity, efficiency and accuracy. In MLPG, the moving least squares approximation is used to construct the trial functions [9]. The obtained shape function does not have the delta function properties. As such, the boundary conditions are difficult to be enforced. Some special techniques have been used to overcome this difficulty in using MLPG, i.e. the Lagrange multiplier method [9], the penalty method [14], and the constraint MLS method [18]. However, the Lagrange multiplier method leads to the unbanded nonpositive definite stiffness matrix. This increases the difficulty to solve the discrete equations. The penalty method requires a proper choice of the penalty factor, which is difficult for some practical problems to choose. The constraint MLS method is very time-consuming in computation. At present, the direct point interpolation seems to be the most efficient method to implement the essential boundary condition in MLPG. The LPIM is the method that employs a set of points to present the problems. The polynomial interpolation functions are constructed with the delta function property. Then, the essential boundary conditions can be easily imposed. In LPIM, a weak form is implemented locally by the weighted residual method, based on the idea of MLPG [12]. It is a truly meshless method, in which the interpolation does not use any element and the integral does not require any background mesh. Regarding to the choice of the test function (weight function), MLPG and LPIM often use the weight function constructed through the same algorithm as the construction of the trial function, based on the Galerkin formulation. In fact, the higher order spline function can also be employed as the weight function. According to LPIM concept, a variation of local point interpolation method (v LPIM) is presented here, in which the weight function is implemented by high-order spline functions, instead of the Galerkin formulation. Using the present v LPIM with the Timoshenko beam theory, the meshless analysis of MEMS devices is conducted for deformations due to the applied voltage between the beam-like component and bottom electrode in the MEMS device. The considered MEMS is a nonlinear system, in which the external force imposed on the system depends on not only the applied voltage but also the deflection of the beam-like component. When the applied voltage reaches a certain value, called the pull-in voltage, the beam-like component will have a contact with the bottom electrode. Both static and dynamic responses of MEMS devices are analyzed in this paper. In the dynamic analysis, the Newmark technique is used to solve the nonlinear dynamic equations. Several MEMS devices are studied and the simulation results are compared with the experiment data and the numerical results by other computational techniques, which demonstrates the efficiency and accuracy of the present v LPIM.
2. Meshless point interpolation technique For the development of the meshless method, the trial function is constructed by the interpolation or approximation. The polynomial point interpolation technique developed by Liu and Gu [10] is used in this paper. One of the advantages of this technique is that it is able to construct the shape function with the delta function properties, and easily imposes essential boundary conditions. Consider a function wðxÞ defined in domain V discretized by a set of field nodes. The interpolation function wðxÞ from the surrounding nodes of a point x using the polynomial basis can be written as wðxÞ ¼
n X
pi ðxÞai ¼ pT ðxÞa
ð1Þ
i¼1
where pi ðxÞ is a monomial in the space coordinates xT ¼ {x y}; n is the number of nodes in the interpolation domain for a interpolation point x; and ai is the coefficient for pi ðxÞ corresponding to the interpolation point x: In the matrix form, we have a ¼ { a1
a2
pT ðxÞ ¼ { 1
··· x
x2
an }T ···
ð2Þ xn21 }
ð3Þ
The coefficients ai in Eq. (1) can be determined by enforcing Eq. (1) at n nodes surrounding the point x: It can be written in the following matrix form w e ¼ P0 a
ð4Þ
where we ¼ { w1
w2
w3
···
wn }T
P0 ¼ { pðx1 Þ pðx2 Þ pðx3 Þ · · · 8 9 > > 1 x1 x21 · · · xn21 1 > > > > > > > > 2 n21 > < 1 x2 x2 · · · x2 > = ¼ . . . . . > .. > .. .. .. .. > > > > > > > > > > : 2 n21 ; 1 xn xn · · · xn
ð5Þ pðxn Þ }
ð6Þ
From Eq. (4), the coefficients are obtained a ¼ P21 0 we
ð7Þ
Then, we have wðxÞ ¼ FðxÞwe
ð8Þ
where FðxÞ is the shape function defined by FðxÞ ¼ pT ðxÞP21 0 ¼ { f1 ðxÞ f2 ðxÞ f3 ðxÞ · · ·
fn ðxÞ }
ð9Þ
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
1263
It should be noted that shape functions possess the delta function property, i.e.
fi ðx ¼ xi Þ ¼ 1;
i ¼ 1; 2; …; n
ð10Þ
fj ðx ¼ xi Þ ¼ 0;
j–i
ð11Þ
n X
fi ðxÞ ¼ 1
ð12Þ Fig. 1. The microswitch simplified as a fixed–fixed beam.
i¼1
The main difficulty in the polynomial PIM is the possible singularity of P0 . Several methods have been proposed to overcome this problem [6,10,19]. However, for the 1-dimensional (1D) problems, there is usually no P0 singularity. In this paper, the MEMS devices are simplified as 1D beam structures.
3. Meshless formulation of the MEMS devices
The boundary conditions are as follows 8 > wðx0 Þ ¼ w; on Gw > > > > > > uðx0 Þ ¼ u; on Gu > > > <
›u
on GM Mðx0 Þ ¼ EI ¼ M; > > › x x¼x0 > > > >
> > > Vðx Þ ¼ GAk u þ ›w
on GV > ¼ V; : 0 s ›x x¼x0
ð15Þ
3.1. Governing equation of the MEMS devices In the simulation of the MEMS devices, electrostatic pull-in characteristic is a well-known sharp instability in the behavior of an elastically supported structure subjected to parallel-plate electrostatic actuation [20]. In this paper, the static and dynamic behavior of the microswitches and the microoptoelectromechanical systems (MOEMS) device are studied subject to the electrostatic voltage. These devices can be simplified as 1D beam structures. Based on Ref. [20] and the theory of elasticity [21], the governing equation for the Timoshenko beam is written as 8
›2 w › ›w > > > < rA ›t2 2 ›x GAks ›x þ u 2f ¼ 0 ; in domain V
> ›2 u › ›u
> ›w > : rI 2 2 EI þu ¼0 þGAks ›x ›x ›x ›t ð13Þ
where Gw ; Gu ; GM and GV are the boundaries of the deflection w; rotation u; moment M; and force V prescribed, respectively. For different boundary constraints, Eq. (15) has different forms. For example, for the fixed end of the beam at x0 ; it becomes (
wðx0 Þ ¼ 0; on Gw
uðx0 Þ ¼ 0; on Gu
ð16Þ
For the free end of the beam at x0 ; Eq. (15) becomes 8 > > <
›u
Mðx0 Þ ¼ EI: x¼x0 ¼ 0;on GM ›x
›w > > : Vðx0 Þ ¼GAks u þ ¼ 0;on GV
›x x¼x0
ð17Þ
The initial conditions at the initial time t0 can be written as
where w is the beam deflection, u the rotation, r the mass density, E the modulus of elasticity, I the moment of inertia, A the cross-section area, G the shear modulus, and ks the shear correction coefficient. f is the nonlinear force obtained from Ref. [20], and it takes the following form f ¼2
10 V~ 2 w g ~ 1þ0:65 w ~ 2g2
ð14Þ
where 10 is the permittivity of vacuum, V~ applied voltage, w ~ beam width, and g the gap between the beam and bottom electrode, g ¼ g0 2wðx;tÞ; where g0 is the distance between the beam and electrode before the beam deflects, as shown in Fig. 1.
8 wðx;t0 Þ ¼ w0 ðxÞ; > > > > > uðx;t0 Þ ¼ u0 ðxÞ; > > < ›wðx;t0 Þ ¼ w^ 0 ðxÞ; in domain V > > › t > > > > > : ›uðx;t0 Þ ¼ u^ ðxÞ; 0 ›t
ð18Þ
3.2. Meshless formula of the MEMS devices The local weak form of Eq. (13), covering a local support domain Vs bounded by Gs ; can be obtained by using
1264
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
the weighted residual technique ) 8ð (
›2 w › ›w > > GAks þ u 2 f dV ¼ 0; v rA 2 2 > > < Vs ›x ›x ›t
#
> ð " ›2 u > › ›u ›w > > EI þ u dV ¼ 0 v rI 2 2 þ GAks : ›x ›x ›x ›t Vs ð19Þ where v is the weight function. Integrating Eq. (19) by parts, we have
function is adopted. It can be written as the following form 8 2 3 4 > < 126 di þ8 di 23 di ; 0 # d # r i v rv rv rv vi ðxÞ ¼ > : 0; di $ r v
where di ¼ lxi 2x0 l is the distance from node xi to point x0 ; and rv is the size of the local support domain. This weight function has higher order continuity. And it can be easily constructed with zero value on the boundary of
8
ð dv ð >
›2 w ›w ›w >
¼0 > GAk þ þ v r A d V þ u V 2 n vGAk u 2 vf d > s s
2 < Vs ›x ›x ›t Vs dx Gs
> 2 ð ð > › u dv ›u ›w ›u
> > þ u dV 2 n v EI v rI 2 dV þ ¼0 EI þ vGAks : ›x ›x ›x Gs ›t Vs Vs dx where n is the unit outward normal to the domain Vs : The boundary Gs for the support domain Vs is usually composed of five parts [16]: the internal boundary Gsi ; the boundaries Gsw ; Gsu ; GsM ; and GsV ; over which the essential boundary conditions w; u and natural boundary conditions M; V are specified. The boundaries Gsw with GsV and Gsu with GsM are mutually disjoint. In Eq. (20), two types of functions, the trial functions and the weight function, are required in the present method. The trial function can be obtained through Eq. (8). The weight function also plays an important role in the implement of the present method. In different meshless methods, there are various forms of the weight function [6]. In the present v LPIM method, the fourth-order spline
ð21Þ
ð20Þ
the local support domain, and thus decreases the computational cost. 3.3. Discritization equations of the MEMS devices In the Timoshenko beam theory, the deflection w and rotation u are independent variables, so they are separately discretized in the space domain, that is (
wðx; tÞ ¼ Fw ðxÞwe ðtÞ
uðx; tÞ ¼ Fu ðxÞue ðtÞ
ð22Þ
where Fw ðxÞ and Fu ðxÞ are shape functions corresponding to the deflection and the rotation, respectively. They can be
Fig. 2. Static deflection of the fixed–fixed beam under different applied voltages.
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
1265
The elements of matrices M; K; and vector f can be written as ð rAvi Fwj dV m11 ij ¼ m12 ij
¼
m22 ij ¼ kij11 ¼ kij12
ð Vs
ð
¼0
rIvi Fuj dV
dFwj dvi dFwj
dV þ n GAks vi GAks dx dx dx Gsi þGsw þGsu þGsM Vs
dvi u u
¼ GAks F dV þ ½nGAks vi Fj dx j Gsi þGsw þGsu þGsM Vs ð
kij21 ¼
kij22
Vs m21 ij
ð Vs
GAks vi
dFwj dV dx
! dvi dFuj u þ GAks vi Fj dV ¼ EI dx dx Vs " # dFuj
þ n EIvi
dx ð
Gsi þGsw þGsu þGsV
fiw ¼
ð Vs
G vi f dV þ ½nvi Vl sV
G fiu ¼ ½nvi M sM The nonlinear dynamic system Eq. (23) is solved by the Newmark method, in which it is assumed that the acceleration varies linearly within time interval ðt; t þ DtÞ: The formulation is given as [22] u_ tþDt ¼ u_ t þ ½ð1 2 dÞu€ t þ du€ tþDt Dt utþDt ¼ ut þ u_ Dt þ
Fig. 3. Dynamic response of the fixed–fixed beam. (a) V~ ¼ 4:0 V; (b) V~ ¼ 8:0 V; (c) V~ ¼ 12:0 V:
1 2 b u€ t þ bu€ tþDt Dt2 2
ð24Þ ð25Þ
The above Newmark method is an implicit method, and it is unconditionally stable provided
d $ 0:5 and b $ 0:25ðd þ 0:5Þ
ð26Þ
computed from Eq. (9). we ðtÞ and ue ðtÞ are the nodal deflections and rotations at time t: Using Eq. (22), the discretized form of Eq. (20) can be given € þ KuðtÞ ¼ fðtÞ MuðtÞ
ð23Þ
where M and K are the mass matrix and stiffness matrix, respectively. uðtÞ is the vector of nodal deflections and € rotations, i.e. uðtÞ ¼ {w1 u1 · · · wn un }T : uðtÞ is the accelerator vector, and fðtÞ is the vector of external force. In this dynamic system, the damping influence is neglected.
Fig. 4. The pull-in process of the fixed–fixed beam under different applied voltages.
1266
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
Fig. 5. Relationship of the critical pull-in voltage and the initial gap for the fixed–fixed beam.
For the simplification, d ¼ 0:5 and b ¼ 0:25 are used in this paper. Because Eq. (23) is a nonlinear equation, an iteration technique is required in each time step. The iteration criteria is vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uX un t ðwiþ1 ð27Þ 2 wij Þ2 # e j
Two microswitches are simplified as the fixed –fixed beam and the cantilever beam, respectively. The MOEMS device is composed of two suspensions connected by a mirror. However, the connections between the suspensions and the mirror are insular. In this device simulation, the suspensions are simplified as cantilever beams, and the mirror is simulated by a bar.
j¼1
4.1. The fixed – fixed beam where n is the number of nodes used, wij and wjiþ1 are the deflections of the ith and ði þ 1Þth iteration steps, respectively. e is a specified accuracy tolerance. However, the nonlinear problem (Eq. (23)) considered in this paper is different from the general nonlinear problems. In the general nonlinear problems, the material properties that are usually nonlinear lead to the mass matrix M or stiffness matrix K be not constant. In this paper, the mass matrix M and stiffness matrix K are constant, but the force acting on the structure is the function of the deflection, i.e. nonlinear loading. Therefore, Eq. (23) is nonlinear and the iteration technique is required in each time step.
The static and dynamic behaviors of the fixed – fixed beam as shown in Fig. 1 are studied firstly. The parameters of this beam are taken as [23]: 80 mm long, 10 mm wide and 0.5 mm thick. The initial gap, g0 ; between the beam and the bottom electrode is 0.7 mm. The Young’s modulus E is 169 GPa, the Poisson ratio is 0.3, the shear correction coefficient is 0.833, and the mass density r is 2231 kg/m3. In analyzing this fixed –fixed beam, 41 distributed nodes are used. In Eq. (13), when the time items are omitted, the static governing equation of the beam can be obtained. Fig. 2 plots the result of the static deflection along the beam under different applied voltages. As the applied voltage increases, the deflection of the beam increases, and the gap
4. Results and discussions Two microswitches and one MOEMS device are simulated by the present meshless v LPIM method.
Fig. 6. The microswitch simplified as a cantilever beam.
Fig. 7. Static deflection of the cantilever beam under different applied voltages.
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
Fig. 8. Dynamic response of the cantilever beam.
Fig. 9. The pull-in process of the cantilever beam under different applied voltages.
between the beam and the bottom electrode decreases. When the applied voltage increases to one certain value, the beam becomes unstable and the center of the beam just touches the bottom electrode. This process is defined
1267
as the pull-in behavior and the ‘certain value’ is defined as the quasi-static critical pull-in voltage. In the implementation of the program, it is defined that if the deflection of the center of the beam, called peak deflection, equals to the initial gap, the corresponding critical pull-in voltage is obtained. In this example, the critical pull-in voltage is 15.1 V. Compared with the critical pull-in value 15.17 V obtained through the experiment and simulation in Ref. [23], the error between these two values is only 0.5%. In the dynamic simulation for the fixed –fixed beam, the time step is taken as 1 £ 1023 ms. Fig. 3 shows the dynamic response of the center of the beam, i.e. peak deflection, under different applied voltages. From this figure, two main remarks can be concluded: one is that when the applied voltage increases, the peak deflection of the beam increases nonlinearly; the other one is that when the applied voltage increases, the fundamental frequency of the beam decreases. Similarly, in the dynamic simulation, one dynamic critical pull-in voltage can be defined. Through calculation, the dynamic critical pull-in voltage of the fixed – fixed beam is 13.8 V. This value matches the result in Ref. [23] for the same beam. Furthermore, from the static and dynamic results, we can find that the quasi-static critical pull-in voltage is larger than the dynamic critical pull-in voltage by about 9%. Fig. 4 plots the dynamic pull-in process at different applied voltages. It can be found that as the applied voltage exceeds the dynamic critical pull-in voltage, the beam collapses onto the bottom electrode very quickly. Fig. 5 demonstrates the relationship between the critical pull-in voltage and the initial gap for both static and dynamic analyses of the beam. It is shown that with the increase of the initial gap, both quasi-static critical pull-in voltage and the dynamic critical pull-in voltage increase. When the initial gap is less than 0.5 mm, the difference between the quasi-static and the dynamic critical pull-in voltages is very small. However, as the initial gap exceeds 0.5 mm, this difference becomes larger.
Fig. 10. Relationship of the critical pull-in voltage and the initial gap for the cantilever beam.
1268
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
Fig. 11. The MOEMS device. (a) Before deflection; (b) after deflection.
4.2. The cantilever beam The considered cantilever beam is shown in Fig. 6. The left end of the beam is fixed and the right end is free. The dimensional and the material parameters for the cantilever beam are identical to the fixed – fixed beam in Section 4.1. And 41 nodes are employed to discretize the cantilever beam. Fig. 7 shows the static deflection along the beam for different applied voltages. As the applied voltage increases, the deflection of the beam increases. Furthermore, the deflection of the free end increases largely. When the applied voltage reaches a certain value, the free end of the beam touches the bottom electrode. Similarly, this voltage is defined as the quasi-static critical pull-in voltage. This value for the considered cantilever beam is 2.33 V. The dynamic simulation is also completed for the cantilever beam. The time step in this computation is taken as 5 £ 1023 ms. The peak deflection, which is at the right end of the cantilever beam, under different applied voltages is studied, as shown in Fig. 8. From the figure, it can be found that the peak deflection of the beam increases nonlinearly as the applied voltage increases. The fundamental frequency of the cantilever beam decreases with the increase of the applied voltages. The dynamic critical pull-in voltage for the cantilever beam is 2.12 V. The error between this value and the quasi-static critical pull-in value is about 9%. The dynamic pull-in
process at different applied voltages is shown in Fig. 9. The same conclusion as the fixed – fixed beam is made, that is, as the applied voltage exceeds the dynamic critical pull-in voltage, the beam collapses onto the bottom electrode very quickly. The relationship between the critical pull-in voltage and the initial gap for both the static and dynamic analyses of the cantilever beam is shown in Fig. 10. The static and dynamic critical pull-in voltages increase as the initial gap increases. When the initial gap is less than 0.5 mm, the difference between the quasi-static and dynamic critical pull-in voltages is very small. However, as the initial gap increases, this difference
Fig. 12. Relationship of the tilt angle and the applied voltage for the MOEMS device.
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270 Table 1 The results of the peak deflection, tilt angle and the mirror strain under different applied voltages for the MOEMS device Applied voltage (V)
Deflection (mm)
Tilt angle (8)
Mirror strain
1.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00
0.000013 0.001304 0.005274 0.012080 0.022055 0.035772 0.054264 0.079578 0.116931 0.204788
0.00149 0.14954 0.60434 1.38402 2.52569 4.09221 6.19402 9.04319 13.16273 22.27279
0.3383 £ 1029 0.3406 £ 1025 0.5563 £ 1024 0.2918 £ 1023 0.9724 £ 1023 0.2556 £ 1023 0.5872 £ 1022 0.1259 £ 1021 0.2698 £ 1021 0.8063 £ 1021
becomes larger. This conclusion is also similar to that for the fixed – fixed beam.
1269
5. Conclusions A variation of local point interpolation method (v LPIM) has been developed for meshless analysis of microelectromechanical devices. It employs higher continuity spline functions as weight function, instead of the Galerkin formulation of the local point interpolation method, and uses the weighted residual technique to discretize the devices. Several MEMS devices are simulated by the Timoshenko beam theory to demonstrate the efficiency of the presently developed meshless method. The numerical results for the nonlinear static and dynamic studies of the MEMS devices are compared with those obtained through experiments and other simulation techniques. It is observed the present meshless v LPIM method in analyzing the MEMS devices is accurate, easy to implement and very efficient.
4.3. The MOEMS device References Another example is the MOEMS device shown in Fig. 11. In this device, two suspensions simplified as cantilever beams are connected by a mirror. The mirror is simplified as a bar and is assumed to be linear elastic. When the voltage is applied, the two beams will deflect in opposite directions, and the free ends of the beam are subjected to external forces that are induced by the bar extension. Following the deflections of two suspensions, the tilt angle of the mirror will be changed. This can be used to change the direction of the reflected light. The parameters employed in this device are: each suspension 10 mm long, 1 mm wide and 0.5 mm thick. The Young’s modulus is 169 GPa, the Poisson ratio is 0.3, and the initial gap between the suspension and the bottom electrode or up electrode is 0.5 mm. The mirror is 1 mm long, and its width, thickness and the Young’s modulus are the same as the suspension. Due to the symmetry of the device, only one suspension and half of the mirror connecting with this suspension are analyzed. In the simulation by the present v LPIM method, the beam is discretized by 41 nodes. Because the beam is analyzed in 1D model, so it only deforms in the vertical direction, and no deformation in the transversal direction (along the axial line). However, the mirror will extend and produce strains along its axial line. This strain will induce a nonlinear force that imposes on the free end of the beam. Fig. 12 presents the curve of the tilt angle and the applied voltage. It can be seen that the tilt angle of the mirror increases with the increase of the applied voltage, and the increase is nonlinear. Table 1 lists the peak deflection, the tilt angle, and the mirror strain at different voltages. From the table, we can conclude that the demand for the different tilt angle of the mirror can be satisfied through changing the applied voltage.
[1] Lyshevski SE. MEMS and NEMS systems, devices, and structures. Boca Raton, FL: CRC Press; 2002. [2] Senturia SD. CAD challenges for microsensors, microactuators, and microsystems. Proc IEEE 1998;86:1611–26. [3] Hung ES, Senturia SD. Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs. IEEE J Microelectromech Syst 1999;8(3):280 –9. [4] Ng TY, Li H, Cheng JQ, Lam KY. A new hybrid meshless-differential order reduction (hM-DOR) method with applications to shape control of smart structures via distributed sensors/actuators. Engng Struct 2003;25:141 –54. [5] Liu GR, Dai KY, Lim KM, Gu YT. A radial point interpolation method for simulation of two-dimensional piezoelectric structures. Smart Mater Struct 2003;12:171–80. [6] Liu GR. Meshfree methods: moving beyond the finite element method. Boca Raton, FL: CRC Press; 2002. [7] Gingold RA, Monaghan JJ. Smooth particle hydrodynamics: theory and applications to non-spherical stars. Monthly Notices R Astronom Soc 1977;181:375–89. [8] Liszka T. An interpolation method for an irregular net of nodes. Int J Numer Meth Engng 1984;20:1599–612. [9] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Meth Engng 1994;37:229–56. [10] Liu GR, Gu YT. A point interpolation method for two-dimensional solid. Int J Numer Meth Engng 2001;50:937–51. [11] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an over-view and recent developments. Computer Meth Appl Mech Engng 1996;139:3–47. [12] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22: 117–27. [13] Atluri SN, Cho JY, Lim HG. Analysis of thin beams, using the meshless local Petrov–Galerkin (MLPG) method, with generalized moving least squares interpolation. Comput Mech 1999;24:334– 47. [14] Gu YT, Liu GR. A meshless local Petrov–Galerkin (MLPG) formulation for static and free vibration analysis of thin plates. CMESComputer Modeling in Engng & Sciences 2001;2(4):463–76. [15] Gu YT, Liu GR. A meshless local Petrov–Galerkin (MLPG) method for free and forced vibration analyses for solids. Comput Mech 2001; 27(3):188– 98.
1270
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
[16] Gu YT, Liu GR. A local point interpolation method for static and dynamic analysis of thin beams. Computer Meth Appl Mech Engng 2001;190:5515–28. [17] Liu GR, Gu YT. A local radial point interpolation method (LR-PIM) for free vibration analyses of 2-D solids. J Sound Vib 2001;246(1):29–46. [18] Liu GR, Yan L. A modified meshless local Petrov–Galerkin method for solids mechanics. Proceeding of Advances in Computational Engineering and Science, Los Angels; 2000. p. 1374–9. [19] Liu GR, Gu YT. A matrix triangularization algorithm for point interpolation method. Proceeding of the Asia-pacific Vibration Conference, Hangzhou China; 2001. p. 1151–4.
[20] Osterberg PM, Senturia SD. M-TEST: a test chip for MEMS material property measurement using electrostatically actuated test structures. J Microelectromech Syst 1997;6(2):107–18. [21] Timoshenko SP, Goodier JN. Theory of elasticity, 3rd ed. New York: McGraw-Hill; 1970. [22] Reddy JN. An introduction to the finite element method, 2nd ed. New York: McGraw Hill; 1993. [23] Ananthasuresh GK, Gupta RK, Senturia SD. An approach to macromodeling of MEMS for nonlinear dynamic simulation. Microelectromech Syst (MEMS), ASME Dyn Syst Control (DSC) Ser 1996;59:401– 7.
A variation of local point interpolation method (v LPIM) for analysis of microelectromechanical systems (MEMS) device Hua Lia,*, Q.X. Wanga, K.Y. Lama,b a
Institute of High Performance Computing, National University of Singapore, 1 Science Park Road, No. 01-01 The Capricorn, Singapore Science Park II, Singapore, Singapore 117528 b Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore, Singapore 119260 Received 28 May 2003; revised 18 August 2003; accepted 18 August 2003 Available online 16 March 2004
Abstract In this paper, a variation of local point interpolation method (v LPIM) is presented for meshless analysis of microelectromechanical systems (MEMS) devices, in which the weight function is assigned by high-order spline functions, instead of the Galerkin formulation of the local point interpolation method (LPIM). The present v LPIM with Timoshenko beam theory are then used to analyze MEMS. In the conventional MEMS design and simulation, time-consuming mesh generation is required. Recently, a number of meshless techniques have been applied to MEMS simulations since they do not require mesh generation. As one of them, the LPIM is based on the local weak form of the partial differential equations and easily implements the boundary conditions due to the shape functions constructed by the PIM with the delta function properties. Several MEMS devices are studied by the present v LPIM. The simulated results are compared with those by other simulation approaches, and it is shown that the present meshless method is very efficient and accurate for the analysis of MEMS devices. q 2004 Elsevier Ltd. All rights reserved. Keywords: Microelectromechanical system; Meshless method; Point interpolation; Static and dynamic analyses
1. Introduction In the past years, attention to microelectromechanical systems (MEMS) has a phenomenal growth because of its wide application in automotive systems (transducers and accelerometers), avionics and aerospace (microscale actuators and sensors), manufacturing and fabrication (microsmart robots), medicine and bioengineering (DNA and genetic code analysis and synthesis, drug delivery, diagnostics and imaging), etc. [1]. However, the development of the increasingly complicated MEMS demands sophisticated simulation techniques [2,3]. The traditional finite element methods or computer-aided design (CAD) systems require the generation of the threedimensional mesh, since the MEMS devices typically involve multiple coupled energy domains and media. It is obvious that mesh generation is computationally expensive and mesh refinement is very difficult, especially for the MEMS devices with complicated geometries and * Corresponding author. Tel.: þ 65-64191-561; fax: þ 65-64191-380. E-mail address: [email protected] (H. Li). 0955-7997/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2003.08.008
multiphysic domains. Therefore, it is necessary to find alternatives to the traditional methods. With the development and the application of the meshless approaches, a lot of meshless modelings and simulations for MEMS (such as microswitches, microsensors/microactuators, etc.) [4,5] have been carried out. Generally, the meshless methods are grouped into two categories according to using or not using integration [6]. The first category is the meshless methods, which do not require integration, based on the strong form of the partial differential equations (PDEs). It includes the smooth particle hydrodynamics [7], the generalized finite difference method [8], etc. The second category is the meshless methods based on the weak form of PDEs, including the element-free Garlerkin method [9], the point interpolation method (PIM) [10], and other methods referred to Refs. [6,11]. Furthermore, in order to avoid the global integration, the local weak forms are also used for the development of local meshless techniques, such as the meshless Petrov – Galerkin method (MLPG) [12 – 15] and the family of local point interpolation method (LPIM) [6,16,17].
1262
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
In the application of the local meshless methods, the choice of the trial function and test function is very important to ensure the computational simplicity, efficiency and accuracy. In MLPG, the moving least squares approximation is used to construct the trial functions [9]. The obtained shape function does not have the delta function properties. As such, the boundary conditions are difficult to be enforced. Some special techniques have been used to overcome this difficulty in using MLPG, i.e. the Lagrange multiplier method [9], the penalty method [14], and the constraint MLS method [18]. However, the Lagrange multiplier method leads to the unbanded nonpositive definite stiffness matrix. This increases the difficulty to solve the discrete equations. The penalty method requires a proper choice of the penalty factor, which is difficult for some practical problems to choose. The constraint MLS method is very time-consuming in computation. At present, the direct point interpolation seems to be the most efficient method to implement the essential boundary condition in MLPG. The LPIM is the method that employs a set of points to present the problems. The polynomial interpolation functions are constructed with the delta function property. Then, the essential boundary conditions can be easily imposed. In LPIM, a weak form is implemented locally by the weighted residual method, based on the idea of MLPG [12]. It is a truly meshless method, in which the interpolation does not use any element and the integral does not require any background mesh. Regarding to the choice of the test function (weight function), MLPG and LPIM often use the weight function constructed through the same algorithm as the construction of the trial function, based on the Galerkin formulation. In fact, the higher order spline function can also be employed as the weight function. According to LPIM concept, a variation of local point interpolation method (v LPIM) is presented here, in which the weight function is implemented by high-order spline functions, instead of the Galerkin formulation. Using the present v LPIM with the Timoshenko beam theory, the meshless analysis of MEMS devices is conducted for deformations due to the applied voltage between the beam-like component and bottom electrode in the MEMS device. The considered MEMS is a nonlinear system, in which the external force imposed on the system depends on not only the applied voltage but also the deflection of the beam-like component. When the applied voltage reaches a certain value, called the pull-in voltage, the beam-like component will have a contact with the bottom electrode. Both static and dynamic responses of MEMS devices are analyzed in this paper. In the dynamic analysis, the Newmark technique is used to solve the nonlinear dynamic equations. Several MEMS devices are studied and the simulation results are compared with the experiment data and the numerical results by other computational techniques, which demonstrates the efficiency and accuracy of the present v LPIM.
2. Meshless point interpolation technique For the development of the meshless method, the trial function is constructed by the interpolation or approximation. The polynomial point interpolation technique developed by Liu and Gu [10] is used in this paper. One of the advantages of this technique is that it is able to construct the shape function with the delta function properties, and easily imposes essential boundary conditions. Consider a function wðxÞ defined in domain V discretized by a set of field nodes. The interpolation function wðxÞ from the surrounding nodes of a point x using the polynomial basis can be written as wðxÞ ¼
n X
pi ðxÞai ¼ pT ðxÞa
ð1Þ
i¼1
where pi ðxÞ is a monomial in the space coordinates xT ¼ {x y}; n is the number of nodes in the interpolation domain for a interpolation point x; and ai is the coefficient for pi ðxÞ corresponding to the interpolation point x: In the matrix form, we have a ¼ { a1
a2
pT ðxÞ ¼ { 1
··· x
x2
an }T ···
ð2Þ xn21 }
ð3Þ
The coefficients ai in Eq. (1) can be determined by enforcing Eq. (1) at n nodes surrounding the point x: It can be written in the following matrix form w e ¼ P0 a
ð4Þ
where we ¼ { w1
w2
w3
···
wn }T
P0 ¼ { pðx1 Þ pðx2 Þ pðx3 Þ · · · 8 9 > > 1 x1 x21 · · · xn21 1 > > > > > > > > 2 n21 > < 1 x2 x2 · · · x2 > = ¼ . . . . . > .. > .. .. .. .. > > > > > > > > > > : 2 n21 ; 1 xn xn · · · xn
ð5Þ pðxn Þ }
ð6Þ
From Eq. (4), the coefficients are obtained a ¼ P21 0 we
ð7Þ
Then, we have wðxÞ ¼ FðxÞwe
ð8Þ
where FðxÞ is the shape function defined by FðxÞ ¼ pT ðxÞP21 0 ¼ { f1 ðxÞ f2 ðxÞ f3 ðxÞ · · ·
fn ðxÞ }
ð9Þ
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
1263
It should be noted that shape functions possess the delta function property, i.e.
fi ðx ¼ xi Þ ¼ 1;
i ¼ 1; 2; …; n
ð10Þ
fj ðx ¼ xi Þ ¼ 0;
j–i
ð11Þ
n X
fi ðxÞ ¼ 1
ð12Þ Fig. 1. The microswitch simplified as a fixed–fixed beam.
i¼1
The main difficulty in the polynomial PIM is the possible singularity of P0 . Several methods have been proposed to overcome this problem [6,10,19]. However, for the 1-dimensional (1D) problems, there is usually no P0 singularity. In this paper, the MEMS devices are simplified as 1D beam structures.
3. Meshless formulation of the MEMS devices
The boundary conditions are as follows 8 > wðx0 Þ ¼ w; on Gw > > > > > > uðx0 Þ ¼ u; on Gu > > > <
›u
on GM Mðx0 Þ ¼ EI ¼ M; > > › x x¼x0 > > > >
> > > Vðx Þ ¼ GAk u þ ›w
on GV > ¼ V; : 0 s ›x x¼x0
ð15Þ
3.1. Governing equation of the MEMS devices In the simulation of the MEMS devices, electrostatic pull-in characteristic is a well-known sharp instability in the behavior of an elastically supported structure subjected to parallel-plate electrostatic actuation [20]. In this paper, the static and dynamic behavior of the microswitches and the microoptoelectromechanical systems (MOEMS) device are studied subject to the electrostatic voltage. These devices can be simplified as 1D beam structures. Based on Ref. [20] and the theory of elasticity [21], the governing equation for the Timoshenko beam is written as 8
›2 w › ›w > > > < rA ›t2 2 ›x GAks ›x þ u 2f ¼ 0 ; in domain V
> ›2 u › ›u
> ›w > : rI 2 2 EI þu ¼0 þGAks ›x ›x ›x ›t ð13Þ
where Gw ; Gu ; GM and GV are the boundaries of the deflection w; rotation u; moment M; and force V prescribed, respectively. For different boundary constraints, Eq. (15) has different forms. For example, for the fixed end of the beam at x0 ; it becomes (
wðx0 Þ ¼ 0; on Gw
uðx0 Þ ¼ 0; on Gu
ð16Þ
For the free end of the beam at x0 ; Eq. (15) becomes 8 > > <
›u
Mðx0 Þ ¼ EI: x¼x0 ¼ 0;on GM ›x
›w > > : Vðx0 Þ ¼GAks u þ ¼ 0;on GV
›x x¼x0
ð17Þ
The initial conditions at the initial time t0 can be written as
where w is the beam deflection, u the rotation, r the mass density, E the modulus of elasticity, I the moment of inertia, A the cross-section area, G the shear modulus, and ks the shear correction coefficient. f is the nonlinear force obtained from Ref. [20], and it takes the following form f ¼2
10 V~ 2 w g ~ 1þ0:65 w ~ 2g2
ð14Þ
where 10 is the permittivity of vacuum, V~ applied voltage, w ~ beam width, and g the gap between the beam and bottom electrode, g ¼ g0 2wðx;tÞ; where g0 is the distance between the beam and electrode before the beam deflects, as shown in Fig. 1.
8 wðx;t0 Þ ¼ w0 ðxÞ; > > > > > uðx;t0 Þ ¼ u0 ðxÞ; > > < ›wðx;t0 Þ ¼ w^ 0 ðxÞ; in domain V > > › t > > > > > : ›uðx;t0 Þ ¼ u^ ðxÞ; 0 ›t
ð18Þ
3.2. Meshless formula of the MEMS devices The local weak form of Eq. (13), covering a local support domain Vs bounded by Gs ; can be obtained by using
1264
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
the weighted residual technique ) 8ð (
›2 w › ›w > > GAks þ u 2 f dV ¼ 0; v rA 2 2 > > < Vs ›x ›x ›t
#
> ð " ›2 u > › ›u ›w > > EI þ u dV ¼ 0 v rI 2 2 þ GAks : ›x ›x ›x ›t Vs ð19Þ where v is the weight function. Integrating Eq. (19) by parts, we have
function is adopted. It can be written as the following form 8 2 3 4 > < 126 di þ8 di 23 di ; 0 # d # r i v rv rv rv vi ðxÞ ¼ > : 0; di $ r v
where di ¼ lxi 2x0 l is the distance from node xi to point x0 ; and rv is the size of the local support domain. This weight function has higher order continuity. And it can be easily constructed with zero value on the boundary of
8
ð dv ð >
›2 w ›w ›w >
¼0 > GAk þ þ v r A d V þ u V 2 n vGAk u 2 vf d > s s
2 < Vs ›x ›x ›t Vs dx Gs
> 2 ð ð > › u dv ›u ›w ›u
> > þ u dV 2 n v EI v rI 2 dV þ ¼0 EI þ vGAks : ›x ›x ›x Gs ›t Vs Vs dx where n is the unit outward normal to the domain Vs : The boundary Gs for the support domain Vs is usually composed of five parts [16]: the internal boundary Gsi ; the boundaries Gsw ; Gsu ; GsM ; and GsV ; over which the essential boundary conditions w; u and natural boundary conditions M; V are specified. The boundaries Gsw with GsV and Gsu with GsM are mutually disjoint. In Eq. (20), two types of functions, the trial functions and the weight function, are required in the present method. The trial function can be obtained through Eq. (8). The weight function also plays an important role in the implement of the present method. In different meshless methods, there are various forms of the weight function [6]. In the present v LPIM method, the fourth-order spline
ð21Þ
ð20Þ
the local support domain, and thus decreases the computational cost. 3.3. Discritization equations of the MEMS devices In the Timoshenko beam theory, the deflection w and rotation u are independent variables, so they are separately discretized in the space domain, that is (
wðx; tÞ ¼ Fw ðxÞwe ðtÞ
uðx; tÞ ¼ Fu ðxÞue ðtÞ
ð22Þ
where Fw ðxÞ and Fu ðxÞ are shape functions corresponding to the deflection and the rotation, respectively. They can be
Fig. 2. Static deflection of the fixed–fixed beam under different applied voltages.
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
1265
The elements of matrices M; K; and vector f can be written as ð rAvi Fwj dV m11 ij ¼ m12 ij
¼
m22 ij ¼ kij11 ¼ kij12
ð Vs
ð
¼0
rIvi Fuj dV
dFwj dvi dFwj
dV þ n GAks vi GAks dx dx dx Gsi þGsw þGsu þGsM Vs
dvi u u
¼ GAks F dV þ ½nGAks vi Fj dx j Gsi þGsw þGsu þGsM Vs ð
kij21 ¼
kij22
Vs m21 ij
ð Vs
GAks vi
dFwj dV dx
! dvi dFuj u þ GAks vi Fj dV ¼ EI dx dx Vs " # dFuj
þ n EIvi
dx ð
Gsi þGsw þGsu þGsV
fiw ¼
ð Vs
G vi f dV þ ½nvi Vl sV
G fiu ¼ ½nvi M sM The nonlinear dynamic system Eq. (23) is solved by the Newmark method, in which it is assumed that the acceleration varies linearly within time interval ðt; t þ DtÞ: The formulation is given as [22] u_ tþDt ¼ u_ t þ ½ð1 2 dÞu€ t þ du€ tþDt Dt utþDt ¼ ut þ u_ Dt þ
Fig. 3. Dynamic response of the fixed–fixed beam. (a) V~ ¼ 4:0 V; (b) V~ ¼ 8:0 V; (c) V~ ¼ 12:0 V:
1 2 b u€ t þ bu€ tþDt Dt2 2
ð24Þ ð25Þ
The above Newmark method is an implicit method, and it is unconditionally stable provided
d $ 0:5 and b $ 0:25ðd þ 0:5Þ
ð26Þ
computed from Eq. (9). we ðtÞ and ue ðtÞ are the nodal deflections and rotations at time t: Using Eq. (22), the discretized form of Eq. (20) can be given € þ KuðtÞ ¼ fðtÞ MuðtÞ
ð23Þ
where M and K are the mass matrix and stiffness matrix, respectively. uðtÞ is the vector of nodal deflections and € rotations, i.e. uðtÞ ¼ {w1 u1 · · · wn un }T : uðtÞ is the accelerator vector, and fðtÞ is the vector of external force. In this dynamic system, the damping influence is neglected.
Fig. 4. The pull-in process of the fixed–fixed beam under different applied voltages.
1266
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
Fig. 5. Relationship of the critical pull-in voltage and the initial gap for the fixed–fixed beam.
For the simplification, d ¼ 0:5 and b ¼ 0:25 are used in this paper. Because Eq. (23) is a nonlinear equation, an iteration technique is required in each time step. The iteration criteria is vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uX un t ðwiþ1 ð27Þ 2 wij Þ2 # e j
Two microswitches are simplified as the fixed –fixed beam and the cantilever beam, respectively. The MOEMS device is composed of two suspensions connected by a mirror. However, the connections between the suspensions and the mirror are insular. In this device simulation, the suspensions are simplified as cantilever beams, and the mirror is simulated by a bar.
j¼1
4.1. The fixed – fixed beam where n is the number of nodes used, wij and wjiþ1 are the deflections of the ith and ði þ 1Þth iteration steps, respectively. e is a specified accuracy tolerance. However, the nonlinear problem (Eq. (23)) considered in this paper is different from the general nonlinear problems. In the general nonlinear problems, the material properties that are usually nonlinear lead to the mass matrix M or stiffness matrix K be not constant. In this paper, the mass matrix M and stiffness matrix K are constant, but the force acting on the structure is the function of the deflection, i.e. nonlinear loading. Therefore, Eq. (23) is nonlinear and the iteration technique is required in each time step.
The static and dynamic behaviors of the fixed – fixed beam as shown in Fig. 1 are studied firstly. The parameters of this beam are taken as [23]: 80 mm long, 10 mm wide and 0.5 mm thick. The initial gap, g0 ; between the beam and the bottom electrode is 0.7 mm. The Young’s modulus E is 169 GPa, the Poisson ratio is 0.3, the shear correction coefficient is 0.833, and the mass density r is 2231 kg/m3. In analyzing this fixed –fixed beam, 41 distributed nodes are used. In Eq. (13), when the time items are omitted, the static governing equation of the beam can be obtained. Fig. 2 plots the result of the static deflection along the beam under different applied voltages. As the applied voltage increases, the deflection of the beam increases, and the gap
4. Results and discussions Two microswitches and one MOEMS device are simulated by the present meshless v LPIM method.
Fig. 6. The microswitch simplified as a cantilever beam.
Fig. 7. Static deflection of the cantilever beam under different applied voltages.
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
Fig. 8. Dynamic response of the cantilever beam.
Fig. 9. The pull-in process of the cantilever beam under different applied voltages.
between the beam and the bottom electrode decreases. When the applied voltage increases to one certain value, the beam becomes unstable and the center of the beam just touches the bottom electrode. This process is defined
1267
as the pull-in behavior and the ‘certain value’ is defined as the quasi-static critical pull-in voltage. In the implementation of the program, it is defined that if the deflection of the center of the beam, called peak deflection, equals to the initial gap, the corresponding critical pull-in voltage is obtained. In this example, the critical pull-in voltage is 15.1 V. Compared with the critical pull-in value 15.17 V obtained through the experiment and simulation in Ref. [23], the error between these two values is only 0.5%. In the dynamic simulation for the fixed –fixed beam, the time step is taken as 1 £ 1023 ms. Fig. 3 shows the dynamic response of the center of the beam, i.e. peak deflection, under different applied voltages. From this figure, two main remarks can be concluded: one is that when the applied voltage increases, the peak deflection of the beam increases nonlinearly; the other one is that when the applied voltage increases, the fundamental frequency of the beam decreases. Similarly, in the dynamic simulation, one dynamic critical pull-in voltage can be defined. Through calculation, the dynamic critical pull-in voltage of the fixed – fixed beam is 13.8 V. This value matches the result in Ref. [23] for the same beam. Furthermore, from the static and dynamic results, we can find that the quasi-static critical pull-in voltage is larger than the dynamic critical pull-in voltage by about 9%. Fig. 4 plots the dynamic pull-in process at different applied voltages. It can be found that as the applied voltage exceeds the dynamic critical pull-in voltage, the beam collapses onto the bottom electrode very quickly. Fig. 5 demonstrates the relationship between the critical pull-in voltage and the initial gap for both static and dynamic analyses of the beam. It is shown that with the increase of the initial gap, both quasi-static critical pull-in voltage and the dynamic critical pull-in voltage increase. When the initial gap is less than 0.5 mm, the difference between the quasi-static and the dynamic critical pull-in voltages is very small. However, as the initial gap exceeds 0.5 mm, this difference becomes larger.
Fig. 10. Relationship of the critical pull-in voltage and the initial gap for the cantilever beam.
1268
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
Fig. 11. The MOEMS device. (a) Before deflection; (b) after deflection.
4.2. The cantilever beam The considered cantilever beam is shown in Fig. 6. The left end of the beam is fixed and the right end is free. The dimensional and the material parameters for the cantilever beam are identical to the fixed – fixed beam in Section 4.1. And 41 nodes are employed to discretize the cantilever beam. Fig. 7 shows the static deflection along the beam for different applied voltages. As the applied voltage increases, the deflection of the beam increases. Furthermore, the deflection of the free end increases largely. When the applied voltage reaches a certain value, the free end of the beam touches the bottom electrode. Similarly, this voltage is defined as the quasi-static critical pull-in voltage. This value for the considered cantilever beam is 2.33 V. The dynamic simulation is also completed for the cantilever beam. The time step in this computation is taken as 5 £ 1023 ms. The peak deflection, which is at the right end of the cantilever beam, under different applied voltages is studied, as shown in Fig. 8. From the figure, it can be found that the peak deflection of the beam increases nonlinearly as the applied voltage increases. The fundamental frequency of the cantilever beam decreases with the increase of the applied voltages. The dynamic critical pull-in voltage for the cantilever beam is 2.12 V. The error between this value and the quasi-static critical pull-in value is about 9%. The dynamic pull-in
process at different applied voltages is shown in Fig. 9. The same conclusion as the fixed – fixed beam is made, that is, as the applied voltage exceeds the dynamic critical pull-in voltage, the beam collapses onto the bottom electrode very quickly. The relationship between the critical pull-in voltage and the initial gap for both the static and dynamic analyses of the cantilever beam is shown in Fig. 10. The static and dynamic critical pull-in voltages increase as the initial gap increases. When the initial gap is less than 0.5 mm, the difference between the quasi-static and dynamic critical pull-in voltages is very small. However, as the initial gap increases, this difference
Fig. 12. Relationship of the tilt angle and the applied voltage for the MOEMS device.
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270 Table 1 The results of the peak deflection, tilt angle and the mirror strain under different applied voltages for the MOEMS device Applied voltage (V)
Deflection (mm)
Tilt angle (8)
Mirror strain
1.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00
0.000013 0.001304 0.005274 0.012080 0.022055 0.035772 0.054264 0.079578 0.116931 0.204788
0.00149 0.14954 0.60434 1.38402 2.52569 4.09221 6.19402 9.04319 13.16273 22.27279
0.3383 £ 1029 0.3406 £ 1025 0.5563 £ 1024 0.2918 £ 1023 0.9724 £ 1023 0.2556 £ 1023 0.5872 £ 1022 0.1259 £ 1021 0.2698 £ 1021 0.8063 £ 1021
becomes larger. This conclusion is also similar to that for the fixed – fixed beam.
1269
5. Conclusions A variation of local point interpolation method (v LPIM) has been developed for meshless analysis of microelectromechanical devices. It employs higher continuity spline functions as weight function, instead of the Galerkin formulation of the local point interpolation method, and uses the weighted residual technique to discretize the devices. Several MEMS devices are simulated by the Timoshenko beam theory to demonstrate the efficiency of the presently developed meshless method. The numerical results for the nonlinear static and dynamic studies of the MEMS devices are compared with those obtained through experiments and other simulation techniques. It is observed the present meshless v LPIM method in analyzing the MEMS devices is accurate, easy to implement and very efficient.
4.3. The MOEMS device References Another example is the MOEMS device shown in Fig. 11. In this device, two suspensions simplified as cantilever beams are connected by a mirror. The mirror is simplified as a bar and is assumed to be linear elastic. When the voltage is applied, the two beams will deflect in opposite directions, and the free ends of the beam are subjected to external forces that are induced by the bar extension. Following the deflections of two suspensions, the tilt angle of the mirror will be changed. This can be used to change the direction of the reflected light. The parameters employed in this device are: each suspension 10 mm long, 1 mm wide and 0.5 mm thick. The Young’s modulus is 169 GPa, the Poisson ratio is 0.3, and the initial gap between the suspension and the bottom electrode or up electrode is 0.5 mm. The mirror is 1 mm long, and its width, thickness and the Young’s modulus are the same as the suspension. Due to the symmetry of the device, only one suspension and half of the mirror connecting with this suspension are analyzed. In the simulation by the present v LPIM method, the beam is discretized by 41 nodes. Because the beam is analyzed in 1D model, so it only deforms in the vertical direction, and no deformation in the transversal direction (along the axial line). However, the mirror will extend and produce strains along its axial line. This strain will induce a nonlinear force that imposes on the free end of the beam. Fig. 12 presents the curve of the tilt angle and the applied voltage. It can be seen that the tilt angle of the mirror increases with the increase of the applied voltage, and the increase is nonlinear. Table 1 lists the peak deflection, the tilt angle, and the mirror strain at different voltages. From the table, we can conclude that the demand for the different tilt angle of the mirror can be satisfied through changing the applied voltage.
[1] Lyshevski SE. MEMS and NEMS systems, devices, and structures. Boca Raton, FL: CRC Press; 2002. [2] Senturia SD. CAD challenges for microsensors, microactuators, and microsystems. Proc IEEE 1998;86:1611–26. [3] Hung ES, Senturia SD. Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs. IEEE J Microelectromech Syst 1999;8(3):280 –9. [4] Ng TY, Li H, Cheng JQ, Lam KY. A new hybrid meshless-differential order reduction (hM-DOR) method with applications to shape control of smart structures via distributed sensors/actuators. Engng Struct 2003;25:141 –54. [5] Liu GR, Dai KY, Lim KM, Gu YT. A radial point interpolation method for simulation of two-dimensional piezoelectric structures. Smart Mater Struct 2003;12:171–80. [6] Liu GR. Meshfree methods: moving beyond the finite element method. Boca Raton, FL: CRC Press; 2002. [7] Gingold RA, Monaghan JJ. Smooth particle hydrodynamics: theory and applications to non-spherical stars. Monthly Notices R Astronom Soc 1977;181:375–89. [8] Liszka T. An interpolation method for an irregular net of nodes. Int J Numer Meth Engng 1984;20:1599–612. [9] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Meth Engng 1994;37:229–56. [10] Liu GR, Gu YT. A point interpolation method for two-dimensional solid. Int J Numer Meth Engng 2001;50:937–51. [11] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an over-view and recent developments. Computer Meth Appl Mech Engng 1996;139:3–47. [12] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22: 117–27. [13] Atluri SN, Cho JY, Lim HG. Analysis of thin beams, using the meshless local Petrov–Galerkin (MLPG) method, with generalized moving least squares interpolation. Comput Mech 1999;24:334– 47. [14] Gu YT, Liu GR. A meshless local Petrov–Galerkin (MLPG) formulation for static and free vibration analysis of thin plates. CMESComputer Modeling in Engng & Sciences 2001;2(4):463–76. [15] Gu YT, Liu GR. A meshless local Petrov–Galerkin (MLPG) method for free and forced vibration analyses for solids. Comput Mech 2001; 27(3):188– 98.
1270
H. Li et al. / Engineering Analysis with Boundary Elements 28 (2004) 1261–1270
[16] Gu YT, Liu GR. A local point interpolation method for static and dynamic analysis of thin beams. Computer Meth Appl Mech Engng 2001;190:5515–28. [17] Liu GR, Gu YT. A local radial point interpolation method (LR-PIM) for free vibration analyses of 2-D solids. J Sound Vib 2001;246(1):29–46. [18] Liu GR, Yan L. A modified meshless local Petrov–Galerkin method for solids mechanics. Proceeding of Advances in Computational Engineering and Science, Los Angels; 2000. p. 1374–9. [19] Liu GR, Gu YT. A matrix triangularization algorithm for point interpolation method. Proceeding of the Asia-pacific Vibration Conference, Hangzhou China; 2001. p. 1151–4.
[20] Osterberg PM, Senturia SD. M-TEST: a test chip for MEMS material property measurement using electrostatically actuated test structures. J Microelectromech Syst 1997;6(2):107–18. [21] Timoshenko SP, Goodier JN. Theory of elasticity, 3rd ed. New York: McGraw-Hill; 1970. [22] Reddy JN. An introduction to the finite element method, 2nd ed. New York: McGraw Hill; 1993. [23] Ananthasuresh GK, Gupta RK, Senturia SD. An approach to macromodeling of MEMS for nonlinear dynamic simulation. Microelectromech Syst (MEMS), ASME Dyn Syst Control (DSC) Ser 1996;59:401– 7.