Ruth V. Sabariego received the M.Sc. degree in Telecommunication Engineering from the University of Vigo, Spain, in 1998. She is currently a Ph.D. student at the Department of Electrical Engineering and Computer Science, University of Liège, Belgium. Her research interest include numerical techniques for computational electromagnetics, applied mathematics, radiation hazards and bio-electromagnetics. Johan Gyselinck graduated in Electrical and Mechanical Engineering at the Ghent University in 1991 and received the Ph.D. degree from the same university in 2000. Since 2000 he is a Postdoctoral Researcher at the Department of Electrical Engineering and Computer Science, University of Liège, Belgium. His research mainly deals with computational electromagnetics, in particular the finite element analysis of electrical machines. Patrick Dular received the electrical engineer degree and Ph.D. degree in applied science from the University of Lige, Belgium, in 1990 and 1994, respectively. Here is currently with the Department of Electrical Engineering and Computer Science, University of Lige, as a Research Associate with the Belgian National Fund for Scientific Research (F.N.R.S.). His research mainly deals with the numerical modeling of coupled problems in electromagnetic systems, including both physical models and numerical techniques coupling. Jeroen De Coster was born in Leuven, Belgium in 1978. He received his MS degree in electromechanical engineering from the Katholieke Universiteit Leuven in 2001. In the same year, he joined the ESAT-MICAS laboratory at the K.U.Leuven as a research assistant. He is currently working towards a PhD, focusing on design, electrical and mechanical characterization of (RF) MEMS devices. Francois Henrotte graduated in engineering (Mechanics-Physics) in 1991 and received the PhD in Applied Sciences in 2000, at the University of Liege, Belgium. In 2000-2004, he was a postdoctoral researcher with the K.U.Leuven, Belgium. Since April 2004, he is with the Institut of Electrical Machines in Aachen, Germany. His main research interests are the numerical modelling of electromechanical systems, the theoretical analysis of electromagnetic forces, differential geometry and applied mathematics. Kay Hameyer received his MS degree in electrical engineering in 1986 from the University of Hanover, Germany. He received his PhD degree from the University of Technology, Berlin, Germany, in 1992. From 1986 to 1988 he worked with the Robert Bosch GmbH in Stuttgart, Germany. In 1988 he became a 1 member of staff at the University of Technology, Berlin, Germany. Beginning of 1993 he was a scientific consultant working on several industrial projects. From 1993 to March 1994 he held a HCM-CEAM fellowship financed by the European Community at the Katholieke Universiteit Leuven, Belgium. From April 1994 to December 2003, he was professor with the KU Leuven. Since January 2004, he is professor with the Institut of Electrical Machines in Aachen, Germany. His research interests are numerical field computation, the design of electrical machines, in particular permanent magnet excited machines and numerical optimisation strategies. 2 Coupled mechanical-electrostatic FE-BE analysis with FMM acceleration–Application to a shunt capacitive MEMS switch R. V. Sabariego, J. Gyselinck, P. Dular, Dept. of Electrical Engineering and Computer Science (ELAP), University of Liège Sart Tilman Campus, Building B28, B-4000 Liège, Belgium J. De Coster, Dept. of Electrical Engineering (ESAT-MICAS), Katholieke Universiteit Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium F. Henrotte, K. Hameyer Dept. of Electrical Engineering (ESAT-ELECTA), Katholieke Universiteit Leuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium Keywords Fast multipole method, Finite element method, Boundary element method, Hybrid method, Microelectromechanical systems Abstract This paper deals with the coupled mechanical-electrostatic analysis of a shunt capacitive MEMS switch. The mechanical and electrostatic parts of the problem are modelled by the FE method and BE method, respectively. The fast multipole method is applied to reduce the storage requirements and the computational cost of the BE electrostatic model. An adaptive truncation expansion of the 3D Laplace Green function is employed. The strong interaction between the mechanical and electrostatic systems is taken into account iteratively. 1 Introduction Electrostatic parallel-plate actuators are widely used in many types of microelectromechanical systems (MEMS). MEMS switches can be used in series or shunt mode and their contacts can be resistive or capacitive (Brown 1998, Tilmans 2002). A shunt capacitive MEMS switch consists of a metal armature (bridge) suspended over a bottom conductor, e.g. the center conductor of a coplanar waveguide, mechanically anchored and electrically connected to the ? The research was carried out in the frame of the Inter-University Attraction Poles IAP P5/34 for fundamental research funded by the Belgian government. P. Dular is a Research Associate with the Belgian National Fund for Scientific Research (F.N.R.S.). ∗ Corresponding author. Email address: R.Sabariego@ulg.ac.be (R. V. Sabariego). Preprint submitted to Elsevier Science 18 February 2004 ground. A thin dielectric film is deposited on the bottom conductor (see Fig. 1). When the bridge is up, the capacitance of the switch is very small and the RF signal freely passes through (the RF switch is on). By applying a bias voltage the switch is actuated: an electrostatic force occurs between top and bottom conductors and the bridge is pulled down, the capacitance increases and causes an RF short to ground (the RF switch is off) (Brown 1998, Tilmans 2002).  !"#9;:<( *'+=,.'>?02324"5      
78    78
77   6   6   78 
78    6   6 

              !"#!$%'&)( *'+-,.'+/+102324"5      
78    78 
78    6   6 

        
8    78 8 
78 Figure 1. Electrostatically actuated capacitive shunt switch implemented on a CPW transmission line. Side and top views These actuators can be treated, in first approximation, as lumped spring-mass systems with a single mechanical degree of freedom (Tilmans 2002). This analysis is helpful for physical insight, but disregards important effects such as the bending of the top plate and the stiction between the bridge and the bottom contact (Brown 1998). The performance of RF MEMS switches strongly depends on the deformation of the top electrode. A detailed knowledge of the exact deformation for an accurate estimate of the capacitance is thus crucial. A boundary element (BE) approach is particularly suited for the analysis of the real electrostatic problem (Farina 2001). Indeed, the BE method provides a rigorous treatment for open problems and allows to consider the deformation without any remeshing. The elastic deformation of the top plate (and the suspension beams) can be handled by means of a finite element (FE) model. It depends directly on the electrostatic force exerted on the bridge and the material properties. The electrostatic field induces a force distribution the value of which increases when the distance between the top and bottom plate diminishes. This interaction between the electrostatic and mechanical systems can be taken into account iteratively. A significant disadvantage of the BE electrostatic model is that it leads to a fully populated system matrix limiting the size of the problems to be handled. The fast multipole method (FMM) (Rokhlin 1983), combined with an iterative solver, e.g. GMRES (Saad 1986), can be employed to overcome this limitation by diminishing both the storage requirements and the computational time. The FMM method has succesfully been applied to solve electrostatic problems in (Nabors 1991, Buchau 2000). In this paper, we discuss the coupled mechanical-electrostratic analysis of a capacitive MEMS shunt switch. Section II outlines the electrostatic BE model 4 of the actuator. The FMM is briefly described in the following section. An adaptive truncation scheme for the 3D Laplace Green function is employed. Section IV deals with the elastic deformation FE model. In the following section, the application example is considered. Simulated results obtained by means of different software packages are briefly compared. 2 Electrostatic BE model We consider an electrostatic problem in 3 . The conductors are embedded in multiple homogeneous isotropic dielectrics and set to fixed potentials. The surfaces of conductors and dielectrics Γ = ΓC ∪ ΓD are discretised with plane triangles. The surface charge density q is assumed to be piecewise constant. The conductors can be replaced by their charge density on their surfaces qc and the homogeneous dielectrics by the polarisation charge qp . The total charge on the interface conductor-dielectric ΓC is given by the sum of both types of charges. Analogously, on the surface between two dielectrics ΓD the total charge is the sum of the polarisation charge due to both dielectrics (Rao 1984). The following system of nq linear equations has to be solved M Q = B,  where Q = q1 . . . qnq  T (1) contains the charge densities on the elements and T depends on the boundary conditions. For an element on B = b 1 . . . b nq the surface of a conductor ΓC , the entry in B is the imposed potential; for an element on the interface between two dielectrics ΓD , the entry in B is zero. The elements of the dense nonsymmetric matrix M when k is an element on a conductor are given by Mk,l = 1 I G(ρk ) dΓ0 ε 0 Γl with G(ρk ) = 1 . 4πρk (2) G(ρk ) is the 3D Laplace Green function, ρk = |r k − r0 | being the distance between a source point r 0 (on Γl ∈ Γ) and an observation point r k (on ΓC ). Considering the continuity of the normal component of the dielectric displacement d = ε e at the dielectric-to-dielectric interface ΓD , the elements of M if k is an element on ΓD read: Mk,l =                εk2 − εk1 ε0 (εk1 + εk2 ) I Γl gradG(ρk ) · nk dΓ0 , k 6= l , (3) 1 , 2ε0 k = l. where nk is the outward-normal unit vector pointing into the dielectric with 5 permittivity εk2 . The integrals in (2) and (3) can be evaluated analytically (Graglia 1993). The electrostatic force F e distribution can be calculated as F e (r) = 1 q(r) e(r) . 2 (4) The electric field e as r approaches the interface conductor-dielectric can be expressed as (Rao 1984): e± (r) = ± n 1 q(r) + h gradG(ρ), q(r0 )iΓ , 2ε0 4πε0 (5) where + indicates the outer face of the conducting surface and − the inner one, n is the normal unit vector pointing outside the conductor and h·, ·iΓ denotes a surface integral on Γ of the product of its arguments. As inside the conductor e− = 0, considering (5), it follows that e+ on the surface of the conductor is given by q(r) . (6) e+ (r) = n ε0 Substituting (6) in (4), the expression of F e as a function of the charge distribution is obtained 1 2 F e (r) = q (r) n . (7) 2 ε0 3 Fast Multipole Method The implementation of the fast multipole method (FMM) requires the groupS ing of the elements on the surface boundary Γ = #g g=1 Γg . A good choice is a scheme based on cubes, i.e. un octree (Nabors 1991, Buchau 2000). Note that in a single level FMM, as described in the present paper, only the finest level of the octree is taken into account. The interactions between distant groups are then determined by means of the multipole expansion of the Laplace Green function. 3.1 Multipole expansion Let Γs be a source group with center r sc and a source point r s , and Γo an observation group with center r oc and an observation point ro . We define the vectors r = r o − roc = (r, θ, φ), r c = r oc − r sc = (rc , θc , φc ) and r 0 = r sc − rs = (r 0 , θ 0 , φ0 ). Omitting the factor 1/4π, the 3D Laplace Green function (2b), with ρ = |ro − r s |, is expanded as (Rokhlin 1983):  X ∞ m X ∞ X u X 1 Dm,n Tm+u,n+v Au,v , =< ρ m=0 n=−m u=0 v=−u with Dm,n (r) = r m Lm n (θ, −φ) , (m + n)! 6 (8) (9) Tm+u,n+v (r c ) = (m + u − (n + v))! n+v Lm+u (θc , φc ), rcm+u+1 Au,v (r 0 ) = r 0u Lvu (θ 0 , −φ0 ) , (u + v)! (10) (11) n ınφ where Lm , Pmn being the Legendre function of degree n (θ, φ) = Pm (cos θ) e m and order n. The imaginary number is denoted ı and < indicates the real part. In practice, the multipole expansion (8) must be truncated by taking 0 ≤ m ≤ p and 0 ≤ u ≤ p, where the truncation number p must be sufficiently large to limit the error to a prescribed value . In most cases, the conventional choice p = log2 (1/) (Rokhlin 1983) is too conservative. Indeed, if r 0  rc and r  rc , a smaller number of terms suffices. Let us consider the radii of the source and observation groups, Rs = maxΓs (r 0 ), Ro = maxΓo (r), and the distance between their centers d. A more economic law p = p(Rs /d, Ro /d, ), proposed by some of the authors in (Sabariego 2004), takes those distances into account. The function grad G in (3) can be expanded in a similar way. It suffices to derive (9) with respect to the coordinates of the observation point. 3.2 Application to the BE model Two groups Γs and Γo are said to be ‘far’ groups if Rs /d < τ and Ro /d < τ , where d is the distance between the group centers and where τ is chosen smaller than 1/2. For demonstrating the FMM, the BE dense matrix M (2-3) can be formally written as M ≈ M near +M far = M near + #g X #g X far Mo,s . (12) o=1 s=1 | {z } Γo , Γs far Let us consider the degrees of freedom qk and ql of q in the respective far groups Γo ∈ ΓC and Γs∈ Γ. Substituting (8) in (2), the contribution to the far corresponding element Mo,s in M far is given by k,l  < p m X X D Mo,k,m,n m=0 n=−m p X u X T A , Mm+u,n+v Ms,l,u,v u=0 v=−u 7  (13) D Mo,k,m,n = Z A Dm,n dΓ , Ms,l,u,v = Γo,k T Mm+u,n+v = Z Au,v dΓ , (14) Γs,l 1 Tm+u,n+v . 4πε0 (15) The iterative solution of the system of algebraic equations requires the multiplication of M far by a trial vector Q. Group by group, the field produced by the electric charge q in the considered group is aggregated into its center by (14b). This aggregated field is then subsequently translated to the centers of all the far groups by (15), and finally the aggregated and translated field is disaggregated into the degrees of freedom of the far groups thanks to (14a). The multiplication M far Q is further accelerated by means of the adaptive truncation scheme following the law p = p(Rs /d, Ro /d, ) (Sabariego 2004). In case of preconditioning of the iterative solver, the preconditioner is based on the sparse matrix comprising the BE near-field interactions. The assembly stage of the FMM consists in calculating and storing the required D T A complex numbers Mo,k,m,n , Mm+u,n+v and Ms,l,u,v . The matrix M far itself is never built. The integrations in (14) are done numerically but as we are dealing with far interactions a limited number of Gauss integration point suffices. The matrix M near is calculated in the conventional way (see previous section) and stored using a sparse storage scheme. For the M D and M A data of a given group, the truncation number p considered during the FMM assembly stage is determined by its closest far group, p = pmax . For the M T data, the truncation number p is determined by the two groups Γs and Γo involved in the translation, p = pso . During the iterative process, the aggregation step is carried out with p = pmax , while p = pso suffices for the translation and disaggregation. 4 Elastic deformation–FE model The upper electrode is deformed by the electrostatic force exerted on it. The elastic equation has to be considered alongside the electrostatic equations. For linear elastic isotropic materials, it reads: DT E D u + F = 0 , (16) where D is the differential operator matrix with transpose DT , E is the elasticity tensor, u is the displacement vector and F is the total force exerted. The elasticity tensor E relates the stress tensor with the strain tensor. It depends on the Young’s modulus E and the Poisson’s ratio ν (Pilkey 2002). 5 Application example The shunt capacitive MEMS switch represented in Fig. 2 is chosen as test case. It concerns a perforated top plate (thickness = 4µm) suspended by a set of 8 
  
    Figure 2. Geometry of shunt capacitive MEMS switch: L c = 475µm, bc = 275µm, Lin = 485µm, bin = 285µm, Ls = 625µm, Lb = 205µm, bb = 20µm and da = 80µm beams, and a bottom plate (thickness = 0.5µm) coated with a thin dielectric layer (thickness = 0.2µm, εr = 7). The beam suspension allows a vertical movement with respect to the fixed bottom plate. The top plate is perforated to facilitate the under-etching of the structure. The dimension of the holes is 25µm × 25µm, with a pitch of 50µm. The mechanical material constants of the top plate are E = 70 Gpa and ν = 0.3. The BE method with FMM acceleration is applied for solving the electrostatic problem while the mechanical problem is handled by a FE model. All the above mentioned methods are implemented in GetDP (GetDP 2003). The behaviour of the switch is simulated using a discretisation consisting of 6544 triangles and 11151 tetrahedra, which yields 6544 degrees of freedom for the piecewise element constant charge q and 56331 degrees of freedom for the second order interpolation of the displacement u. The optimal number of FMM groups (for this particular mesh) is found to be 35. The maximum and average truncation number are pmax = 6 and pav = 4 for Rf ar = 135µm and  = 10−6 . The electrostatic and mechanical systems are solved iteratively obtaining the new electrostatic force distribution and the new displacement. The number of iterations required for sufficient convergence of e.g. the capacitance increases as the applied bias voltage approaches the pull-in voltage and the deformation of the top plate becomes bigger. The calculated zero-voltage capacitance CV =0 and pull-in voltage VIN are 0.36 pF and 14.2 V, respectively. The deformation of the top electrode for a bias voltage of 11 V for the successive iterations is shown in Fig. 3. Convergence is achieved after 9 iterations. 9 -0.05 -0.1 u (µm) -0.15 -0.2 -0.25 -0.3 -0.35 -0.4 -300 -200 -100 0 100 200 300 x(µm) Figure 3. Convergence of the vertical displacement along a line through the suspension beams and perforated plate for an applied bias voltage of 11 V The results obtained with GetDP (GetDP 2003) are compared with those given by the commercial software packages Coventor (Coventor 2003) and FemLab (Femlab 1997-2004). In the simulations performed with the commercial programs, only a quarter of the geometry is considered. In the Coventor simulation, the electrostatic part is modelled by means of the BE method while the mechanical part is dealt with using the FE method and second order elements. Symmetry boundary conditions are only taken into account for the mechanical problem. In the FemLab computation, the whole electromechanical problem is solved with the FE method. Symmetry conditions are imposed for the electrostatic problem. With regard to the mechanical part, the elastic behaviour of the suspension (beams) is approximated by a stiffness constant (Brown 1998, Tilmans 2002). For the face of the top electrode that is coupled with the suspension, the displacement is obtained by dividing the total electrostatic force by the stiffness constant. The nominal capacitance CV =0 obtained by Coventor and Femlab is 0.4pF and 0.37pF, respectively. The pull-in calculated voltage is 14.24 V for Coventor and 17.25 V. The computed value of the capacitance as a function of the applied voltage is represented in Fig. 4 for the three different solvers. The curves C − V obtained with GetDP and FemLab agree well for low applied voltage, when the deformation is small. As the applied voltage increases, an accurate estimate of the displacement becomes critical, the approximation used for the suspension does not suffice. On the contrary, the agreement between the curves obtained with GetDP and Coventor is better as the voltage increases. The influence of three quarters of the device are disregarded for the electrostatic computation but the mechanical part is solved accurately. Fig. 5 shows the maximum vertical displacement of the top electrode as a function of the applied bias voltage. A good agreement between the values obtained by means of GetDP and Coventor is observed. The use of the approximation for the mechanical problem with a stiffness constant for modelling the suspension in the FemLab 10 0.6 GetDP CoventorWare FemLab Capacitance (pF) 0.55 0.5 0.45 0.4 0.35 0 2 4 6 8 10 12 14 16 18 Voltage (V) Figure 4. Calculated capacitance vs the applied bias voltage simulation explains the divergence of the curves. Maximum displacement (µm) 1.4 GetDP CoventorWare FemLab 1.2 1 0.8 0.6 0.4 0.2 0 0 2 4 6 8 10 12 14 16 18 Voltage (V) Figure 5. Maximum vertical displacement of the top electrode vs the applied bias voltage 6 Conclusion A shunt capacitive MEMS switch has been modelled. The BE method, accelerated by the FMM, and the FE method have been applied to solve the electrostatic and mechanical problem, respectively. An adaptive truncation scheme for the 3D Laplace Green function has been employed. The results have been compared with those obtained with the commercial packages Coventor and FemLab. References Brown, E. R. (1998), “RF-MEMS switches for reconfigurable integrated circuits”, IEEE Trans. on Microwave Theory and Techniques, Vol. 46, No. 11, November, pp. 1868–80. Buchau, A., Huber, C. J., Rieger W. and Rucker W. M. (2000), “Fast BEM computations with the adaptive multilevel fast multipole method”, IEEE Transactions on Magnetics, Vol. 36, No. 4, July, pp. 680–4. 11 Coventor, Inc. (2003), Accelerating MEMS Innovations, 4001 Weston Parkway, Cary, NC 27513, USA, http://www.coventor.com. Farina, M. and Rozzi, T. (2001), “A 3D integral equation-based approach to the analysis of real-life MMICS–Application to Microelectromechanical systems”, IEEE Trans. on Microwave Theory and Techniques, Vol. 49, No. 12, pp. 2235– 40. Femlab, Multiphysics modelling (1997-2004), COMSOL, Inc., 1100 Glendon Avenue, 17th Floor Los Angeles, CA 90024, USA, http://www.femlab.com. GetDP: a General Environment for the Treatment of Discrete Problems (2003), Department of Electrical Engineering (ELAP), Institut Montefiore, University of Liège, Sart Tilman Campus, Building B28, B-4000 Liège, Belgium, http://www.geuz.org/getdp/. Graglia, R. D. (1993), “On the numerical integration of the linear shape functions times the 3D Green’s Function of its gradient on a plane triangle”, IEEE Trans. on Antennas and Propagation, Vol. 41, No. 10, pp. 1448–55. Nabors, K. and White, J. (1991), “FastCap: A multipole accelerated 3-D capacitance extraction program”, IEEE Transactions on Computer-Aided Design, Vol. 10, No. 11, pp. 1447–59. Pilkey, Walter D. (2002), Analysis and design of elastic beams: Computational methods, John Wiley & Sons, New York. Rao, S. M., Sarkar, T. K. and Harrington,R. F. (1984), “The electrostatic field of conducting bodies in multiple dielectric media”, IEEE Trans. on Microwave Theory and Techniques, Vol. 32, No. 11, pp. 1441–8. Rokhlin, V. (1983), “Rapid solution of integral equations of classical potential theory”, Journal of Computational Physics, Vol. 60, pp. 187–207. Saad, Y. and Schultz, M. H. (1986), “GMRES: A Generalized Minimal Residual Algorithm for solving nonsymmetric linear systems”, SIAM J. Sci. Comput., Vol. 7, No. 3, pp. 856–69. Sabariego, R. V., Gyselinck, J., Dular, P., Geuzaine, C. and Legros W. (2004), “Fast Multipole Acceleration of the Hybrid Finite Element–Boundary Element Analysis of 3D Eddy Current Problems”, accepted for publication in IEEE Transactions on Magnetics. Tilmans, H. A. C. (2002),“MEMS Components for Wireless Communications (invited paper)”, Proc. EuroSensors XVI, 15-18 September, Prague, Czech Republic, pp. 1–34 (CD-ROM only). 12