Eugenio Brusa
email: eugenio.brusa@uniud.it

Antonio Della-Gaspera

Mircea Gh

Munteanu
email: 2munteanu@uniud.it

Antonio Della Gaspera

Validation of compact models of microcantilever actuators for RF-MEMS application

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Reconfigurable layout is a typical feature of circuits used for radio-frequency application [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF]. This goal is achieved by means of cantilever and double clamped microbeams, electrostatically actuated, being used as switches, resonators and varactors [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF][START_REF] Bona | Microsystems Mechanical Design[END_REF][START_REF] Senturia | Microsystems design[END_REF][START_REF] Najar | Modeling and design of variable-geometry electrostatic microactuators[END_REF]. Design needs for a precise prediction of pull-in condition and frequency response. An effective modelling of these static and dynamic behaviours is rather difficult, because of the electromechanical coupling. Electromechanical forces nonlinearly depend on mechanical displacement, electric charge and voltage. In presence of large strain or displacement a structural nonlinear solution has to be implemented [START_REF] Choi | Improved analysis of microbeams under mechanical and electrostatic loads[END_REF][START_REF] Pamidighantam | Pull-in voltage analysis of electrostatically actuated beam structures with fixedfixed and fixed-free end conditions[END_REF][START_REF] Somà | Large deflections of microbeams under electrostatic loads[END_REF]. Analytical solutions for microcantilever and double clamped microbeams were formulated and included the effects of stretching and large displacement, i.e. the so-called geometrical nonlinearity [START_REF] Choi | Improved analysis of microbeams under mechanical and electrostatic loads[END_REF][START_REF] Pamidighantam | Pull-in voltage analysis of electrostatically actuated beam structures with fixedfixed and fixed-free end conditions[END_REF]. Corrections suitable to predict fringing effect of electric field were even proposed [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF][START_REF] Senturia | Microsystems design[END_REF]. These models assume an ideal beam geometry, which differs from the actual structure in some details of the constraint and electrode shape. Microcantilevers currently proposed by industry look like specimens depicted in Figs. 1 and2. A first geometry is based on "in-plane" bending test, i.e. microbeam deflection occurs in a plane parallel to the profiling system's target. In "out-ofplane" bending actuation microbeam tip moves towards the target [START_REF] Tsou | On the out-of-plane deformation of Vshaped micromachined beams[END_REF]. Often microbeam is a part of a wider structure, e.g. it connects the rigid plate of a varactor to the fixed frame [START_REF] Gaddi | Low-loss ohmic RF-MEMS switches with interdigitated electrode topology[END_REF]. All these aspects motivate the need for a library of structural numerical models within the simulator used to predict the response of the whole electronic circuit. Accuracy and fastness are main requirements for this modelling activity. Authors demonstrated in previous papers that static pull-in of in-plane bending specimens is accurately predicted by a numerical solution based on sequential approach with voltage increments [START_REF] Somà | Large deflections of microbeams under electrostatic loads[END_REF][START_REF] Munteanu | Geometrical nonlinearities of electrostatically actuated microbeams[END_REF][START_REF] Ballestra | Experimental characterization of electrostatically actuated in-plane bending of microcantilevers[END_REF]. Geometric nonlinearity becomes relevant for large displacements, close to pull-in voltage. A double nonlinear solution, including geometrical effect, has to be implemented. In this case iterative solution is applied at each step of the sequential approach and makes extremely heavy the computational effort. Iteration can be avoided by using a special finite beam element and a non-incremental formulation [START_REF] Munteanu | Geometrical nonlinearities of electrostatically actuated microbeams[END_REF][START_REF] Gh | A new FEM approach to the coupledfield analysis of electrostatic microactuators dynamics[END_REF][START_REF] Brusa | Coupled-field FEM nonlinear dynamics analysis of continuous microsystems by non incremental approach[END_REF]. Effectiveness of 2D FEM models in case of inplane bending specimens was a little bit surprising, because no electric force concentration due to the finite dimensions of the geometry was taken into account [START_REF] Ballestra | Experimental characterization of electrostatically actuated in-plane bending of microcantilevers[END_REF]. In this paper a detailed analysis of in-plane microcantilevers is performed to complete that investigation. Results are compared to those of out-ofplane microcantilevers, which show an initial curvature of microbeam. Main features of the two above mentioned layouts are described. To define suitable criteria to proceed with a model order reduction useful for the dynamic analysis, limits of the coupled field analysis in 2D models are discussed, by investigating some three dimensional effects of the electric field. Numerical results are simultaneously compared to the experiments done on microspecimens made of epitaxial polysilicon and gold. Experimental validation of numerical models developed in previous papers was performed on eight geometries, described in Table I. FEM static analysis was performed in ANSYS code, as sequential solution with mesh morphing in dielectric region, then through the non-incremental FEM sequential approach tested in [START_REF] Munteanu | Geometrical nonlinearities of electrostatically actuated microbeams[END_REF][START_REF] Gh | A new FEM approach to the coupledfield analysis of electrostatic microactuators dynamics[END_REF][START_REF] Brusa | Coupled-field FEM nonlinear dynamics analysis of continuous microsystems by non incremental approach[END_REF] and by a combined sequential FEM/BEM solution [START_REF] Ballestra | Experimental characterization of electrostatically actuated in-plane bending of microcantilevers[END_REF]. Last two methods were implemented in Matlab. An additional comparison included results of the Discrete Geometric Approach (DGA), recently proposed [START_REF] Bettini | Static behaviour prediction of microelectrostatic actuators by discrete geometric approaches (DGA)[END_REF]. As Fig. 3 shows in case of geometry 4, all the 2D approaches converge to the experimental curve, provided that Young modulus of epitaxial polysilicon was clearly identified and geometrical nonlinearity due to the large tip displacement was included. FEM 3D (ANSYS, SOLID122 electrostatic, SOLID185, elastic) analysis revealed that the actual distribution of the electrostatic force is fairly different in terms of average and peak values (Fig. 4). Nevertheless, effectiveness of 2D models appeared surprisingly good. A deeper investigation of three dimensional effects of electric field allowed finding that agreement between 2D models and experiments was assured by constraint and electrode geometries. In case of geometry 5, where half-width is c = 7.5 µm, ratio between the peak values of force computed by 3D FEM model and 2D respectively was 1.88. This result assumed that actuation voltage was applied to the counter-electrode, equal width for electrode and counter-electrode, rounded tip, no surface behind the beam. All these assumptions play a significant role in case of in-plane actuators (Fig. 5). Charge concentration at the tip and along the edges increases electrostatic force, more largely as peak than as average value. This concentration is localized on a small area (Fig. 4). Rectangular and rounded tips have higher force ratio 3D/2D than sharp triangular tip. Indefinitely long counter-electrode and zero voltage applied on the microbeam bring the above force ratio up to 2.1. Actually, wafer surface below the lateral edge of microcantilever decreases this ratio to 1.05, thus allowing 2D model predicting the actual pull-in. In case of equal width for upper and counter-electrode ratio tends to 1, being higher when counter-electrode width is larger than the electrode's one.

Geometric characterization of in-plane microcantilever is rather difficult. Profiling system offers a very high resolution along the optical axis, while it is lower on the target plane. Thickness and gap measurements for in-plane microcantilever is less accurate, as Table I shows. Numerical prediction of pullin is consequently ineffective, if nominal values of these parameters are inputted. For a given correspondence of numerical and experimental values, and in case of geometry 5 following discretizations were set used: DGA (FORTRAN, sequential method) 15000 elements, 32000 DoFs (electrical); 756 elements and 4000 Dofs (structural); FEM (ANSYS; iterative, mesh morphing) 80 PLANE183 (solid beams), 3000 PLANE 121 (electrical); FEM (MATLAB; sequential; non incremental; mesh morphing) 3036 elements, 346 nodes (electrical); 41 nodes, 40 Timoshenko two-node beam elements; FEM/BEM (MATLAB, sequential, non-incremental) 337 two-node boundary elements, 188 nodes (electrical); 31 nodes, 30 two-node Timoshenko beam. As it is clearly described in [START_REF] Ekeom | Electrostatic-Mechanical Direct Coupling and Fem-bem Coupling for CMUT Electroacoustic Parameters Identification[END_REF], BEM allows reducing the number of DoFs in dielectric region, but boundary elements assure a better accuracy within the element field than on boundaries. To predict accurately voltages on the microbeam a mesh refinement is required, although DoFs are less than in FEM.

III. SPECIMENS FOR OUT-OF-PLANE BENDING TEST

Four geometries of golden microcantilevers for out-of-plane bending were built to perform a parametric analysis. Each one included several specimens. Precision in measuring gap and thickness is here higher than in case of in-plane actuators. Nevertheless, specimens exhibit some differences in length, thickness and gap. This layout has two peculiarities. Counter-electrode only partially fills the gap, and the anchor is a structural component with a defined geometry. A crucial aspect was the initial curvature of geometry 9,10 and 11 due to some differences of diffusion of Chromium of the seed-layer among the deposited layers. Models had to include this curvature to fit experimental pull-in voltage. In fact, while microfabrication may induce a residual stress gradient across the beam section in double clamped beam, in microcantilevers stress vanishes since it imposes an initial strain and curvature which moves the free tip. For given initial strain ε 0 , curvature κ 0 , accidental thermal effects and axial or flexural preloads, N 0 and M 0 respectively, stress-strain relations integrated at beam's cross section become:

      +               -             =      
Symbols mean Young modulus, E, beam section, A, beam second moment of area, J, axial effort, N, bending moment, M. Temperature distribution may include a constant contribution T 0 along beam thickness (y axis), and a distribution T(x,y) variable along beam length (x) and thickness (y).

Experimental results were similar for the four geometries tested. Those of geometry 10 are depicted in Fig. 9. FEM and FEM/BEM approaches converge to a numerical solution which overestimates the actual pull-in voltage. In this case benefits of in-plane layout are absent. Fringing effect is more relevant. To fit experiments it was required to perform a 3D FEM analysis and compute the correction factor for the electromechanical force. It was observed that electric field in 2D models, including separately microbeam length and width, gave easily this number. FEM 3D models did show some problems because of mesh morphing operation applied to a so narrow gap. A sensitivity analysis confirmed that electrostatic pressure is quite uniformely distributed along the length and the width of the microbeam, while only peak values of force strongly

Fig. 1 .

 1 Fig. 1. Microcantilever built for in-plane bending test (top view) .

Fig. 2 .

 2 Fig. 2. Microcantilever built for out-of-plane bending test (front view) .

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Example of experimental validation on in-plane microcantilever

9- 11

 11 April 2008 ©EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5

Fig. 6 .

 6 Fig. 6. Potential distribution around the microcantilever half-width in case of null voltage applied to counter-electrode (left) or to beam surface (right).

Fig. 7 .

 7 Fig. 7. 2D and 3D predictions of potential distribution around the microcantilever half-width in case of null voltage applied to beam surface and in presence of wafer surface (right edge in 2D model).

Fig. 8 .

 8 Fig. 8. Effect of finite dimensions of counter-electrode on in-plane microcantilever half-width.

TABLE I MICROCANTILEVER

 I

SPECIMENS FOR IN-PLANE BENDING TEST

  Table II summarizes relevant parameters. Numerical data are written by describing the range of measured values among different specimens and the measurement errors.

TABLE II MICROCANTILEVER

 II SPECIMENS FOR OUT-OF-PLANE BENDING TEST

	n.	Length	Width	Thickness	Gap	y
		[µm]	[µm]	[µm]	[µm]	[µm]
	9	531:535	32:33	2.9:3.0	2.88:2.99	4.15:6.6
		± 0.3	± 0.3	± 0.5.10 -4	± 0.5.10 -4	± 0.5.10 -4
	10 190 ± 0.3	32	1.8	2.97:3.17
			± 0.3	± 0.5.10 -4

ACKNOWLEDGMENT This work was partially funded by the Italian Ministry of University under grant 2005/2005091729, Operative Unit of Udine (p.i. E.Brusa).

©EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5 depend on thickness and gap values. In fact, an electrical analysis on the undeformed configuration of the microbeam allowed computing a correction factor suitable to find experimental pull-in.

IV. DISCUSSION

Model order reduction of nonlinear and second order dynamic microsystems is still a difficult task. No definitive approaches were successfully tested, although some methods demonstrated to be effective in some specific application [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF][START_REF] Bona | Microsystems Mechanical Design[END_REF][START_REF] Senturia | Microsystems design[END_REF][START_REF] Ostergaard | Finite Element based reduced order modelling of micro electro mechanical systems[END_REF][START_REF] Bechtold | Fast simulation of electrothermal MEMS[END_REF].

Nonlinearity is a crucial aspect and microcantilevers exhibit both electromechanical and geometric nonlinearities. Choice is either solving an analytical formulation of the coupled problem, by reducing DoFs involved, or resorting to a numerical sequential solution. In this case mesh morphing inhibits the use of MOR methods [START_REF] Bechtold | Fast simulation of electrothermal MEMS[END_REF]. Ad hoc linearization was already proposed, in absence of geometric nonlinearity [START_REF] Salimbahrami | A simulation free reduction scheme and nonlinear modelling of an electrostatic beam[END_REF], while geometric nonlinear MEMS can be characterized by 2D static models. These allow identifying microsystem stiffness to be used together with damping for dynamic analysis. For in-plane configuration 3D effects of electric field are less relevant, but a very accurate measure of design parameters is required. In case of out-of-plane layout it is just the opposite. Experimental validation demonstrated that force input for 2D models can be calibrated on 3D FEM electrical analyses. Dynamic analyses can be then performed.

V. CONCLUSION

Two dimensional models are often considered poorly effective to predict static and dynamic behaviour of microbeam RF-MEMS. Microcantilevers with in-plane or out-of-plane bending can be accurately and fairly fast analysed by 2D models, based on sequential non incremental approach implemented in FEM, FEM/BEM or DGA. 2D model tuning can be done by performing a FEM analysis of three dimensional effects of electric field. Geometric nonlinearity is relevant for all the specimens tested. In-plane microcantilevers analysis suffers any inaccuracy in measuring the parameters used as inputs. For out-of-plane microbenders fringing and three dimensional effects of electric field have to be carefully evaluated together with the initial curvature, often present.

Where model order reduction techniques fail because of the double nonlinearity of actuation and geometrical effects, 2D models appear suitable to extract few lumped parameters to perform dynamic analysis. This procedure requires an evaluation of electric field singularities to correct the electrostatic force input.


T Veijola

A Lehtovuori

Numerical and Compact Modelling of Squeeze-film Damping in RF MEMS Resonators

Oscillatory gas flow in squeeze-film dampers is studied up to frequencies where the length of the acoustic wave is comparable with the dimensions of the air gap. Damping and spring forces are calculated both numerically and analytically from the linearized 2D Navier-Stokes equations. In addition to the low frequency region of inertialess gas, where the use of the Reynolds equation is limited, the new model considers several additional phenomena. These are the inertia of the gas, the transition from isothermal to adiabatic conditions, and the gap resonances at frequencies where the acoustic wavelength is comparable to the air gap height. Velocity and temperature slip conditions are considered to make the model valid in micromechanical structures where the air gap heights are of the order of a micrometer. An approximate compact model is derived combining the low frequency model and the gap resonance model. The accuracy of the compact model is studied by comparing its response to the numerical results calculated with the finite element method. The agreement is very good in a wide frequency band when the ratio of the damper width and the gap height is greater than 10. The numerical study and the compact model are directly applicable in predicting the damping and the resonance frequency shift due to air in RF MEMS resonators having narrow air gap widths and operating at frequencies where the wavelengths become comparable to the flow channel dimensions.

I. INTRODUCTION

It has been reported RF MEMS resonators that operate at frequencies up to 1 GHz [START_REF] Wang | 1.156-GHz Self-Aligned Vibrating Micromechanical Disk Resonator[END_REF], [START_REF] Wang | 1.51-GHz Nanocrystalline Diamond Micromechanical Disk Resonator With Material-Mismatched Isolating Support[END_REF]. Since these devices can operate also at atmospheric pressures, models predicting the quality factor and resonance frequency shift due to air are needed. The challenges in the modelling are the small acoustic wavelengths that might become comparable with the dimension of the damper width and even the air gap height, and the small dimensions of micromechanical devices that invalidate the continuum flow assumptions.

Squeeze-film damping in air gaps of oscillating structures has been traditionally modelled with the Reynolds equation [START_REF] Blech | On Isothermal Squeeze Films[END_REF], [START_REF] Griffin | A Study of Fluid Squeeze-Film Damping[END_REF]. It considers the viscous gas flow and the compressibility and assumes constant pressure across the air gap. With an effective viscosity, the rarefied gas phenomenon can be considered. This extends the model to be used for micromechanical devices and or low ambient pressure. However, the Reynolds equation is usable only up to a certain frequency where the inertial forces can be neglected (Reynolds number 1). This limitation can be avoided using a modified form of the Reynolds equation considering the inertia in the gap flow [START_REF] Veijola | Compact Models for Squeezed-Film Dampers with Inertial and Rarefied Gas Effects[END_REF] together with non-isothermal conditions [START_REF] Beltman | Air Loads on a Rigid Plate Oscillating Normal to Fixed Surface[END_REF], [START_REF] Beltman | Viscothermal Wave Propagation Including Acousto-Elastic Interaction, Part I: Theory[END_REF]. This viscoelastic wave propagation model is not directly applicable to model MEMS devices, since it assumes continuum boundary conditions.

When the frequency of the oscillation is so high that the length of the acoustic wave becomes comparable to the height of the air gap, the constant pressure assumption across the air gap breaks down. This situation has been studied and modelled in [START_REF] Veijola | Model for Gas Damping in Air Gaps of RF MEMS Resonators[END_REF]. The study shows that since the gas is trapped in the gap, it is justified to assume closed damper borders in deriving the model. The analytic model derived in [START_REF] Veijola | Model for Gas Damping in Air Gaps of RF MEMS Resonators[END_REF] models accurately the damper in the higher frequency region where the air gap resonances occur.

In this paper, numerical simulation of the linearized Navier-Stokes equations are first used to study the oscillating flow in the air gap of a 1D squeeze-film damper in a wide frequency band. Slip velocity and temperature boundary conditions at the surfaces, and zero pressure conditions at the damper borders are applied. A finite element method (FEM) software is used in the study.

Next, reduced linearized Navier-Stokes equations are presented that consider the variable pressure across the air gap to be able to include the gap resonances in the model [START_REF] Veijola | Model for Gas Damping in Air Gaps of RF MEMS Resonators[END_REF]. The gap resonances are solved at high frequencies and the squeeze film phenomenon is solved at lower frequencies. Then, these results are combined to have a compact approximate model for a wide frequency regime. Finally, the compact model is compared with the numerical model to verify its accuracy.

II. OSCILLATING GAS FLOW IN A SQUEEZE-FILM DAMPER

A. Topology of the damper

Figure 1 shows the topology of the studied 1D squeezefilm damper. It is assumed that the 3rd dimension l y is much larger than l x justifying the study of the 2D gas flow in the cross-section. The gas is bounded by two rigid parallel surfaces, and the upper surface is oscillating with a small amplitude of w 0 in the z direction. Trivial zero-pressure boundary conditions are used at the borders of the damper.

B. Damping and spring forces

The damper is characterized with a mechanical impedance Z m = -F w0 , that is calculated from the force F acting on the moving surface having a velocity of w 0 in the z direction. The damping coefficient c is the real part of Z m and the spring coefficient k is the real part of iωZ m .

C. Characteristic numbers

The behaviour of the flow in a narrow air gap is described by a few characteristic numbers. The Reynolds number Re (the square of the shear wave number s) is the ratio between inertial and viscous forces: s 2 = Re = ωh 2 0 ρ 0 /η, where ρ 0 is the density of the gas, η is the viscosity coefficient, and h 0 is the characteristic height of the air gap.

The reduced frequency k = ωh 0 /c 0 is a measure of the ratio between the gap width h 0 and the acoustic wave length. c 0 = γp 0 /ρ 0 is the speed of sound in free space, where γ = c P /c V is the specific heat ratio and p 0 is the ambient pressure. The Prandtl number Pr characterizes the thermal properties. Here we use the square root of Pr,

φ = √ Pr = ηc P κ , (1)
where c P is the specific heat at constant pressure and κ is the thermal conductivity.

A frequency-domain analysis with small perturbation amplitudes is assumed, and the equations are presented in normalized form. The small-amplitude perturbation variables: pressure, temperature and density are specified with p = p 0 + pe iωt , T = T 0 + T e iωt , ρ = ρ 0 + ρe iωt , (2) respectively.

D. Linearized Navier-Stokes equations

Beltman [START_REF] Beltman | Air Loads on a Rigid Plate Oscillating Normal to Fixed Surface[END_REF] has derived linearized dimensionless form of the Navier-Stokes equations. The equations are in 2D (ydirection is not considered):

iu = -g kγ ∂p ∂x + 1 s 2 g 2 ∂ 2 u ∂x 2 + ∂ 2 u ∂z 2 + g 3s 2 ∂ ∂x g ∂u ∂x + ∂w ∂z (3) iw = -1 kγ ∂p ∂z + 1 s 2 g 2 ∂ 2 w ∂x 2 + ∂ 2 w ∂z 2 + 1 3s 2 ∂ ∂z g ∂u ∂x + ∂w ∂z (4) g ∂u ∂x + ∂w ∂z = -ikρ (5) p = ρ + T (6) iT = 1 s 2 φ 2 g 2 ∂ 2 T ∂x 2 + ∂ 2 T ∂z 2 + i γ -1 γ p, ( 7
)
where u and w are velocity components in the x and z directions, respectively. The equations are in normalized form such that velocities u and w are normalized to the speed of sound c 0 and the dimensions x and z are normalized with the characteristic dimensions l x and h 0 , respectively. p, ρ, and T represent small relative variations around p 0 , ρ 0 and T A , respectively. The narrowness of the gap is g = h 0 /l x .

E. Numerical modeling

To study gas flow in the air gap, pressure, velocity, and temperature of the 1D-damper in Fig. 1 are simulated with a FEM solver for visco-acoustic flow [START_REF] Malinen | A Finite Element Method for the Modeling of Thermo-Viscous Effects in Acoustics[END_REF] in multiphysical FEM software ELMER [START_REF] Pironneau | Elmer -Finite Element Solver for Multiphysical Problems[END_REF]. It solves linearized Navier-Stokes equation presented in ( 3)-( 7) as a function of frequency.

In this study, parameters for air at standard atmosphere conditions are used, see Table I. The air gap height is h 0 = 1µm and the length is l x = 5µm. The upper surface oscillates in the z-direction with a constant velocity amplitude of 1 m/s (w 0 = 1). To demonstrate the trapping of the gas in the air gap at high frequencies, Fig. 2 shows forces acting on the upper surface as a function of frequency in two cases. First, the borders are assumed to be open (zero pressure) and then they are assumed to be closed (zero velocity). The amplitude responses are identical above a certain frequency indicating that at high frequencies the flow velocity in the x direction becomes insignificant. Hence, at high frequencies the gas is trapped in the gap corresponding the closed border assumption, even for open border conditions.

Figure 3 shows simulated pressure, velocity, and temperature amplitude distributions for the open border case at certain frequencies, represented with verticals in Fig. 2. These frequencies represent different flow regions that will be studied closer. Resonances occur both in the x and y directions: squeeze film phenomenon at lower frequency when the gas flows out from the gap and gap resonances at higher frequency when the gas is trapped in the gap. (b) First resonance in the x direction. Its frequency f x,N is determined by the wave propagation velocity in the x direction and the length of the structure. The squeeze film resonant frequencies are

f x,N ≈ Nc e 4L , N = 2, 4, . . . ( 8
)
where L is the characteristic length of the damper and c e is the effective speed of sound in the x-direction. The temperature profile in Fig. 3 indicates non-isothermal conditions. (c) First (anti)resonance f z,1 in the x-direction: pressure is not constant across the gap, small velocity (trapped situation). The N th air gap resonance frequency is

f z,N ≈ Nc 0 4h 0 , N = 1, 2, 3, . . . ( 9
)
where h 0 is the height of the gap. Odd values of N give antiresonances, while even values of N give resonances. Air gap resonances occur at clearly higher frequency than the squeeze film resonance if the gap is narrow (h 0 l x ). Non-isothermal conditions. (d) 2nd gap resonance f z,2 : pressure, velocity and temperature vary mainly in the z-direction, practically adiabatic conditions.

F. Reduced equations for the approximative model

In order to derive the compact model, a reduced set of NSequations is needed. The complete set of 2D equations ( 3)-( 7) will be reduced by removing insignificant terms. Equation ( 3) is first written in the following form:

iu = - g kγ ∂p ∂x + 4g 2 3s 2 ∂ 2 u ∂x 2 + 1 s 2 ∂ ∂z ∂u ∂z + g 3 ∂w ∂x ( 10
)
Since velocity u varies less in the x-direction than in the zdirection, and g is small, term ∂ 2 u/∂x 2 can be ignored. Also, due to translational velocity excitation, the term ∂w/∂x is very small and is ignored here. At low frequencies, k is small and (4) reduces simply to ∂p/∂z = 0. At high frequencies ∂ 2 w/∂x 2 is small due to the excitation, and ∂u/∂x is small in (4) due to the trapped situation. Temperature variations are slower in the x direction than in the z direction, and g is small in [START_REF] Beltman | Viscothermal Wave Propagation Including Acousto-Elastic Interaction, Part I: Theory[END_REF], making the term ∂ 2 T /∂x 2 to vanish. The reduced equations are:

iu = - g kγ ∂p ∂x + 1 s 2 ∂ 2 u ∂z 2 (11) iw = - 1 kγ ∂p ∂z + 4 3s 2 ∂ 2 w ∂z 2 (12)
g ∂u ∂x + ∂w ∂z = -ikρ (13) p = ρ + T (14) iT = 1 s 2 φ 2 ∂ 2 T ∂z 2 + i γ -1 γ p. ( 15
)
These equations are identical with the "narrow gap" equations specified in [START_REF] Beltman | Air Loads on a Rigid Plate Oscillating Normal to Fixed Surface[END_REF], but they include an additional equation [START_REF] Veijola | Extending the Validity of Existing Squeezed-Film Damper Models with Elongations of Surface Dimensions[END_REF]. The strategy to have an approximate solution for the equations above is the following. First, the equations are further reduced to have the trapped air situation. This is accomplished assuming flow in the z direction only, that is, u(x) = 0, and solving the velocity profile w(z). Next, the narrow gap situation is assumed ignoring (4), and a differential equation for the pressure distribution is derived.

G. Air gap resonance region

It was shown numerically that at high frequencies the gas is trapped in the air gap, and moves only in the z direction. The problem has been studied closer in [START_REF] Veijola | Model for Gas Damping in Air Gaps of RF MEMS Resonators[END_REF], where both a simplified model with adiabatic conditions and a complete model are given.

H. Squeeze film damping region

Next, the lower frequency (squeeze-film damping) region is studied in order to derive the compact model for the mechanical impedance of the damper. This is accomplished by using the equations ( 11)-( 15) for narrow gap problems, that is, Eq, ( 12) will be ignored in the following analysis.

The solution strategy is similar to the one in [START_REF] Beltman | Air Loads on a Rigid Plate Oscillating Normal to Fixed Surface[END_REF], but here the slip velocity and temperature conditions are used instead of the continuum boundary conditions.

The slip boundary conditions for the gas velocity u and the temperature T are used at the surfaces:

u = ∓K s ∂u ∂z z=0,z=1 T = ∓K T ∂T ∂z z=0,z=1 (16)
The Knudsen number K n = λ/h 0 is a measure of gas rarefaction. It is a ratio between the mean free path λ and the nominal displacement h 0 . Here, the mean free path is specified as

λ = η ρ 0 π 2RT 0 (17)
where R = c Vc P . In this paper, we use the quantity K s = σ P λ/h 0 as a measure of the rarefaction instead of K n . Since λ is inversely proportional to pressure, K s increases as the pressure drops below the ambient pressure. In this paper σ P = 1 is assumed. K T is the "thermal Knudsen number" [START_REF] Karniadakis | Micro Flows, Fundamentals and Simulation[END_REF]

K T = 2 -α T α T 2γ γ + 1 K n φ 2 , ( 18
)
where α T is the energy accommodation coefficient and γ is the specific heat ratio.

Since the pressure was assumed to be constant on a cross section, p and ∂p/∂x are independent of z in ( 11) and (15). Therefore, the temperature can be solved from (15). Solution is searched in form T (z) = b+c cosh[ξ(z - 1 2 )] and applying boundary conditions in ( 16) results in

T = γ -1 γ pA(φs, K T , z), ( 19
)
where the temperature profile is described with

A(s, K s , z) = 1 - cosh[ √ is(z -1 2 )] cosh( √ is/2) + √ isK s sinh( √ is/2) . ( 20
)
Remark, that although the pressure is assumed to be constant in z-direction, the temperature depends on z. Next, the density is solved from ( 14) resulting in ρ = p/m(φs, K T , z) where

m(φs, K T , z) = 1 - γ -1 γ A(φs, K T , z) -1 . ( 21
)
Equation ( 11) has the same form as (15) for temperature and u can be solved when p = p(x) yielding

u = ig kγ ∂p ∂x A(s, K s , z). ( 22
)
The differential equation for the pressure distribution in the x-direction is derived from ( 13) by averaging the velocity u and density ρ across the gap:

g 1 0 ∂u ∂x dz + ik 1 0 ρ dz = - ∂w ∂z z=1 . ( 23
)
Since the mechanical impedance is specified at z = 1, the term in the right hand side has not been averaged across the gap. Instead of the average velocity in the z direction, the excitation velocity at the upper surface (z = 1) is used here. Inserting u from ( 22) and ρ from ( 14) into (23) gives a general frequency-domain model for 1D air gaps

g 2 B(s, K s ) ikγ ∂ 2 p(x) ∂x 2 - ik n(φs, K T ) p(x) = ∂w ∂z z=1 . ( 24
)
where

B(s, K s ) = 1 0 A(s, K s , z) dz (25) = √ is -(2 -is 2 K s ) tanh( √ is/2) √ is + is 2 K s tanh( √ is/2)
and the similar integration of m(φs, K T , z) gives the polytropic constant

n(φs, K T ) = 1 - γ -1 γ B(φs, K T ) -1 . ( 26
)
The function n is the result of thermal effects. It is generally complex-valued signifying a phase shift between the pressure and density perturbation. Hence, it models the losses due to gas compressibility that are important in the transition regime when conditions change from isothermal to adiabatic. At low frequencies (φs 1), the flow is isothermal leading to n = 1 and at high frequencies (φs 1), B(s, K s ) approaches 1 and n → γ (adiabatic conditions).

Equation ( 24) is solved with trivial pressure boundary conditions p(± 1 2 ) = 0. The solution is searched in form where σ c is the complex squeeze number

p(x) = b cosh(cx) + d, yielding p(x) = n(φs, K T ) ik cosh( √ iσ c x) cosh( √ iσ c /2) -1 ∂w ∂z z=1 , (27)
σ c = ik 2 γ g 2 B(s, K s )n(φs, K T ) = is 2 σ 12B(s, K s )n(φs, K T ) ,
(28) where the traditional squeeze number σ is defined as σ = 12ηl 2

x ω/(p 0 h 2 0 ). It specifies the ratio between the spring force due to the gas compressibility and the force due to the viscous flow.

I. Complete approximate model

The resulting force acting on the surface is

F = 1/2 -1/2 1/2 -1/2 p(x) dx dy (29) = n(φs, K T ) ik 2 √ iσ c tanh( iσ c /2) -1 ∂w ∂z z=1 ,
where the velocity excitation is calculated when z = 1.

In (32) the last derivative is still unknown. Two different velocity functions w(z) has been derived in [START_REF] Veijola | Model for Gas Damping in Air Gaps of RF MEMS Resonators[END_REF] assuming closed damper boundaries and gap resonances.

The simpler model gives at z = 1

∂w ∂z z=1 = w 0 q tanh(q) (30)
and the more complicated model gives

∂w ∂z z=1 = C 1 r 1 e r1 + C 2 r 2 e r2 + C 3 r 3 e r3 + C 4 r 4 e r4 , ( 31
)
respectively. Coefficients C i and r i are given in the Appendix.

The un-normalized force becomes

F = - l x l y n(φs, K T )p 0 iωh 0 2 √ iσ c tanh( iσ c /2) -1 ∂w ∂z z=h0 .
(32) Equation ( 32) is the complete model to calculate the force as a function of frequency for a velocity excitation. The damping force is the real part of F .

III. RESULTS

In the following, the results of the compact model are compared to numerical FEM simulations which can be considered as exact results. Then, the importance of the included features in the compact model is studied by comparing the responses of published models in a wide frequency band. In the following figures, magnitude or real part is denoted with × at left-side axis and phase (or spring coefficient) with 2 at right-side axis.

A. 2D FEM Simulations

FEM simulations were performed with a solver for dissipative acoustic flow [START_REF] Malinen | A Finite Element Method for the Modeling of Thermo-Viscous Effects in Acoustics[END_REF] included in Elmer [START_REF] Pironneau | Elmer -Finite Element Solver for Multiphysical Problems[END_REF] software. The solver gives the results directly at specified frequencies. Since the compact model considers a 1D damper, 2D FEM simulations are sufficient here. A sinusoidal velocity amplitude of 1 m/s was used as the excitation. The symmetry of the structure was utilized in the FEM simulations; only a half of the air gap was simulated. Slip velocity and temperature conditions were applied, and the thermal boundary conditions were set for ideally thermally conducting surfaces. Boundary conditions p(±l x /2) = 0 were used at the borders. A mesh of 8000 elements was used, and the simulation was performed at 81 frequencies from 10 kHz to 1 GHz. The grid convergence showed that the element count is sufficient. Gas parameters for air in Table I were used in the simulations.

B. Validity of the compact model

Several numerical simulations were performed to be able to test the validity and accuracy of the compact model for the 1D squeeze-film damper derived in the previous section.

Here, the accurate model for the air gap region in Appendix A is used. Table II summarizes the dimensions of all simulated damper topologies. The simulated topologies were limited to cases where the ratio l x /h 0 ≥ 10. With smaller ratios the zero pressure border condition is not justified. Also, the pressure is no more constant across the gap even at small frequencies. These phenomena have been studied in [START_REF] Veijola | Extending the Validity of Existing Squeezed-Film Damper Models with Elongations of Surface Dimensions[END_REF] regime is between the squeeze-film and the gap resonance regions.

Table II summarizes the resonant frequencies detected from the amplitude/phase responses of all simulated topologies. Here, the resonant frequencies are those where the phase of the impedance crosses zero.

To see the validity of the resonant frequency approximations in [START_REF] Veijola | Model for Gas Damping in Air Gaps of RF MEMS Resonators[END_REF] and [START_REF] Malinen | A Finite Element Method for the Modeling of Thermo-Viscous Effects in Acoustics[END_REF], they are compared with the results in Table II. The effective speed of sound propagating in the gap is

c e = c 0 5(1 + 2K s ) 6γ (33)
The fundamental resonant frequencies f x,2 and f z,1 become 6.94 MHz and 43.52 MHz respectively, for a topology of a width of 20 µm and gap height of 2 µm.

IV. CONCLUSIONS

With these models, the Reynolds equation has been extended to be applicable to rapidly oscillating surfaces and rarefied gas conditions in addition to modelling damping at low frequencies. It was shown that the compact model is in good agreement with numerical FEM simulation results. Slip conditions are used to have an accurate model also for small air gaps typical in MEMS devices.

In this study, idealistic border conditions for the damper borders were assumed. These conditions are not justified in cases where the aspect ratio (l x /h 0 ) is not small. To model accurately such cases, the flow outside the gap should be considered [START_REF] Veijola | Extending the Validity of Existing Squeezed-Film Damper Models with Elongations of Surface Dimensions[END_REF]. It is straightforward to extend the approach to 2D squeeze damper topology considering the y dimension.

The solver for viscous acoustic flow in the Elmer software [START_REF] Pironneau | Elmer -Finite Element Solver for Multiphysical Problems[END_REF] turned out to be very useful in studying the squeeze-film phenomena at the high frequency regime.

Fig. 1 .

 1 Fig. 1. The cross-section of squeeze-film damper. The y dimension of the damper ly is assumed to me much larger than the x dimension lx.

5 9Fig. 3 .Fig. 2 .

 532 Fig. 3. Amplitude profiles of pressure p, velocity u, and temperature T at different frequencies. a) f = 1 MHz, b) f = 30 MHz, c) f = 100 MHz, d) f = 180 MHz

  Finally, the flow profile w(z) is inserted in the differential c EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5 9-11 April 2008 equation that is solved for perpendicular motion and zero pressure boundary conditions.

Figure 4 5 9Fig. 4 .

 454 Figure 4 compares frequency responses of numerical and compact models in a case when l x = 20 µm and h 0 = 2 µm.The simulated data is presented in two forms: (a) amplitude/phase of the impedance Z m , and (b) damping/spring coefficients. Both responses show that the compact model response is very accurate at low and high frequencies, but there is a noticeable error in the compact model just below the first gap resonance frequency. This was expected, since this

APPENDIX A

An exact solution for ( 12) -( 15) is presented considering the non-adiabatic thermal conditions and boundary conditions for velocity and temperature. After some manipulation, the following fourth order differential equation results

(A1) for velocity w and

(A2) for temperature T as well. In this Appendix, the derivatives with the respect to z are denoted with primes (T = ∂ 4 T /∂z 4 ), and the dependence on z is denoted by, e.g., w(z) or T (z). These homogenous linear equations with constant complex coefficients have a characteristic equation

The roots of the characteristic equation are

and thus the solution of (A1) is

Constants C 1 , C 2 , C 3 , and C 4 are determined using boundary conditions:

Therefore, the temperature is written as function of velocity w. Equations ( 12) -( 15) reduce now to

where

)

)

Solving T from (A10) and (A11) and p from ( 5), ( 6) and (A16) yields

where B 1 and B 2 are the auxiliary variables:

Equation (A16) can be used to utilize the boundary conditions for temperature in (A9) to solve the velocity:

After applying boundary conditions for velocity in (A8) in addition to the conditions above, the following system of equations results:

where

Solving the system of equations (A19)-(A22) gives the coefficients C i in (A7) which are

)

where G = w 0 /(e r3e r4 ), (A27)

H 2 = (e r2e r4 )/(e r3e r4 ), (A29)

)

), (A32)


Denis Juillard
email: jerome.juillard@supelec.fr.

Aubry

A Missoffe
email: alexia.missoffe@supelec.fr

J Juillard

D Aubry

A piecewise-linear reduced-order model of squeeze-film damping for deformable structures including large displacement effects

A piecewise-linear reduced-order model of squeeze-film damping for deformable structures including large displacement effects

INTRODUCTION

Correct modelling of damping is essential to capture the dynamic behaviour of a MEMS device. Our interest is squeeze-film damping which models the behaviour of a fluid in small gaps between a fixed surface and a structure moving perpendicular to this surface. The lateral dimensions of the surfaces are large compared to the gap and the system is considered isothermal. Squeeze film damping is then governed by the Reynolds equation [START_REF] Bao | Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes[END_REF]:

3 ( ) 12
G G P P P t µ   ∂ ∇ ∇ =   ∂   (1)
where ( )

t y x G , ,
is the distance between the moving and the fixed surface, ( )

t y x P , ,
is the pressure, µ the effective viscosity of the fluid [START_REF] Bao | Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes[END_REF]. For small excitation frequencies or amplitudes the squeezed film behaves as a linear damper. For larger amplitudes or frequencies, the gas has no time to flow away and the pressure builds up creating a stiffening effect coupled to a nonlinear damping. A complete review on this equation and its different regimes can be found in [START_REF] Bao | Squeeze film air damping in MEMS[END_REF].

Coupling the Reynolds equation to the equation governing the mechanical behaviour of the microstructure leads to a nonlinear system of partial differential equations (PDEs), which has no analytical solution and must be simplified thanks to some assumptions. The most commonly made assumptions are the following:

-uniform displacements, i.e. 0 = ∂ ∂ = ∂ ∂ y G x G -steady-state sinusoidal excitation, i.e. ( ) t G G e ω sin = [4] -small displacements, i.e. g G G + = 0 and 0 G g <<
, where 0 G is the nominal gap of the structure at rest or close to a static equilibrium [START_REF] Younis | Microplate Modeling under Coupled Structural-Fluidic-Electrostatic Forces[END_REF]. -small pressure variations, i.e. p P P + = 0 and 0 P p << , where 0 P is the ambient pressure. These hypotheses prove to be useful in a variety of applications, if only for gaining insight of nonlinear damping phenomena. However, in many cases, it is difficult to justify their use: for example, it is clear to see that none of the first three hypotheses holds when trying to estimate the switching time of a micro-switch. Most micro-switches do not undergo uniform displacements, nor can these displacements be considered small, and the behaviour of a micro-switch is fundamentally transient.

To date, the most notable attempts to tackle the problem of reduced-order modelling (ROM) of squeeze-film damping with large, non-uniform displacements have been made by Younis et al. [START_REF] Younis | Microplate Modeling under Coupled Structural-Fluidic-Electrostatic Forces[END_REF][START_REF] Younis | Simulation of squeeze-film damping of microplates actuated by large electrostatic load[END_REF][START_REF] Younis | A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping[END_REF], Mehner et al. [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF], Yang et al. [START_REF] Yang | Macromodeling of coupled-domain MEMS devices with electrostatic and electrothermal effects[END_REF], Hung and Senturia [START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF] and Rewienski and White [START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]. In [START_REF] Younis | Microplate Modeling under Coupled Structural-Fluidic-Electrostatic Forces[END_REF][START_REF] Younis | Simulation of squeeze-film damping of microplates actuated by large electrostatic load[END_REF][START_REF] Younis | A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping[END_REF], the authors solve the linearized Reynolds equation for a displacement of the beam around an operating point. In [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF], the authors use a modal projection method to calculate modal frequency-dependent damping and stiffening coefficients close to a determined operating point. To extend this approach to large displacements, Mehner [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF] gives an analytical expression of these coefficients as a function of mechanical modal coordinates established by fitting of simulation data for different initial deformations. These approaches are all based on several steady-state sinusoidal calculations [START_REF] Younis | Microplate Modeling under Coupled Structural-Fluidic-Electrostatic Forces[END_REF][START_REF] Younis | Simulation of squeeze-film damping of microplates actuated by large electrostatic load[END_REF][START_REF] Younis | A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping[END_REF] or simulations [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF], which increase the time for setting up the reduced-order model. The most general approaches may well be those developed in [START_REF] Yang | Macromodeling of coupled-domain MEMS devices with electrostatic and electrothermal effects[END_REF][START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF][START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]. These approaches, although very general, have a high computational cost (because of the nonlinear /multiphysics /transient simulation they require) and their accuracy depends, to some degree, on the choice of the training trajectory.

ISBN: 978-2-35500-006-5

                + ∂ ∂ =         ∇ ∇ 0 3 1 12 P p G t p G µ (2)
The first reduction step is based on modal projection of (2) which is first transformed via a change of variable on p (3).

) 2 / 3 ( - = G p ϕ (3)
The aim of this change of variable is to obtain a spatial operator on φ for which the Laplacian eigenmodes are more relevant than for the operator in (2), conserving its selfadjoint property thus guarantying convergence of the solution. The reduced-order model may be written as:

s x H x f s x A ) ( )) ( ) ( ( = - dt d (4) with ∫∫ Ω - Ω = d G P f l l 2 / 1 0 2 ϕ (5)         Ω ∆ + - = ∫∫ Ω d G G P H l k kl k kl ϕ ϕ δ λ µ 2 / 3 2 / 3 2 0 12 (6) and ∫∫ Ω - Ω = d G A l k kl ϕ ϕ 2 (7)
where x and s are respectively the vectors of the mechanical modal coordinates of the moving structure, and the modified squeeze coordinates. For a structure under electrostatic actuation, one may write the full coupled model as: where p e (x) is the electrostatic force. The Euler-Bernoulli equation is used to model the mechanical part. Large displacement effects can be taken into account. Let us write (8) in a more compact way:

( ) ( ) ( ) ( ) ( ) ( )           +                     - =                     -                    
) ( )) ( ( )) ( ( )) ( ( ( t u t t t dt d z B z f z g + = (9)
where z is the state vector including all modal coordinates.

III. PIECEWISE LINEARIZATION BASED ON AN INPUT TRAJECTORY

The cost of evaluating the terms g(z(t)), f(z(t)), and B(z(t)) can be reduced using the piecewise linear approach described in [START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]. The fact that the coefficients in (3) only depend on the mechanical modal coordinates reduces the cost of construction of the piecewise linear model, which has the following structure:

) ( )) ( ) ( )( ( )) ) ( ( ) ( 1 0 1 0 z F z z JF z f z z (z JG z g z i i i elec i s i i i i s i i w w dt d + - + =       - + ∑ ∑ - = - = (10)
where i JG and i JF are respectively the jacobians of the functions g(z) and f(z) at the linearization points i z and (z) F elec corresponds to the electrostatic force. Two problems arise from the piecewise linearization: the choice of the linearization points and the weighting procedure. We choose the linearization points from a simulation of the fully nonlinear reduced model [START_REF] Yang | Macromodeling of coupled-domain MEMS devices with electrostatic and electrothermal effects[END_REF]. A new linearization point is chosen when a point is "far enough" from the already chosen points. The state variables must be normalized to calculate the distances in state-space. The weighting procedure is the one described in [14].We work on the example of a microswitch also treated in [START_REF] Yang | Macromodeling of coupled-domain MEMS devices with electrostatic and electrothermal effects[END_REF][START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF][START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]. We use the piecewise linear reduced order model to determine the switching time of the device for a step voltage between 9 and 10.5 V at atmospheric ambient pressure. The 21 linearization points are chosen along a 9.5V input training trajectory. The time integration scheme is not exactly implicit because the weights are calculated from the last step point. Fig. 1 shows the response to a 10V input using a linear model, the modal projection model, and the piecewise linearized model. Fig. 2 shows the experimental and simulated data presented in [START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF] and results of our piecewise linear reduced model of order 6 for the switching time.

IV. TRAJECTORY-INDEPENDENT REDUCED ORDER MODEL

We can notice that the nonlinear terms in (8) depend only on the mechanical modal coordinate and that the model is linear with respect to the modified squeeze coordinates. It can be written as:

PE(x) z x B x FP z x G + = - ) ( )) ( ) ( ( dt d (11)
The piecewise linearization described in section II doesn't take advantage of this structure. This second piecewise linearized model is based on the linearization of the terms G, FP and B. As they only depend on the on mechanical coordinate, there is no need for a training trajectory. It is sufficient to discretize in an appropriate way the space corresponding to the mechanical coordinates. This is a great advantage as the resulting model does not depend on the relevance of a training trajectory, as opposed to the model presented in section II. On the other hand the resulting equation is nonlinear which increases the resolution cost. This equation is integrated in time using a fully implicit Euler scheme. At each step the nonlinear equation is solved using the Matlab "fsolve" function. Table 1 shows results for an input step voltage of 9.1 V. This model yields better displacement error and correct results are obtained for a minimum number of linearization points of around 15 which is less than the model obtained using the approach described in [START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF].

V. CONCLUSION

We have presented a piecewise linear model of squeezefilm damping for flexible structure and large displacements with the restriction of small pressure variations. The model is based on a modal projection method and is then piecewise linearized using two different methods. The advantages and drawbacks of each method were pointed out. The first method [START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF] is very general whereas the second takes advantage of a specificity of the equation structure. This last method appears to be more efficient, although the final system that must be solved remains nonlinear. The pull-in time of a microswitch was determined and compared to experimental and other simulation data.

  .[START_REF] Bao | Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes[END_REF] has then the following form:

Fig. 1 - 5 Fig. 2 -

 152 Fig. 1-Middle point displacement for the 9.5 V step voltage training trajectory, a 10 V full reduced model simulation, piecewise linearized model, and linear model..

Table 1

 1 Middle point displacement error and pull-in time error for different number of linearization points compared to the full reduced model.

	Number of	Middle point	Pull-in time
	linearization	displacement	error (in %)
	points	error (in %)
	10	21	7
	14	5	2.5
	19	1	2
	30	0.5	0.1


Alberto Ballestra
email: alberto.ballestra@polito.it

Eugenio Brusa
email: eugenio.brusa@uniud.it

Giorgio De Pasquale
email: giorgio.depasquale@polito.it

Mircea Gh

Munteanu
email: 3munteanu@uniud.it

Aurelio Somà
email: aurelio.soma@polito.it

RF-MEMS beam components: FEM modelling and experimental identification of pull-in in presence of residual stress

In this paper an experimental validation of numerical approaches aimed to predict the coupled behaviour of microbeams for out-of-plane bending tests is performed. This work completes a previous investigation concerning in plane microbeams bending.

Often out-of-plane microcantilevers and clamped-clamped microbeams suffer the presence of residual strain and stress, which affect the value of pull-in voltage. In case of microcantilever an accurate modelling includes the effect of the initial curvature due to microfabrication. In double clamped microbeams a preloading applied by tensile stress is considered. Geometrical nonlinearity caused by mechanical coupling between axial and flexural behaviour is detected and modelled.

Experimental results demonstrate a good agreement between FEM approaches proposed and tests. A fairly fast and accurate prediction of pull-in condition is performed, thus numerical models can be used to identify residual stress in microbridges by reverse analysis from the measured value of pull-in voltage.

I.

INTRODUCTION

Microcantilevers and microbridges are currently widely used in RF applications as microswitches and microresonators [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF], [START_REF] Gaddi | Low-loss ohmic RF-MEMS switches with interdigitated electrode topology[END_REF] and in experimental micromechanics, where materials mechanical properties and strength are measured. Therefore, it is required implementing efficient numerical models to predict the electromechanical performance of microstructures actuated by electric field [START_REF] Gaddi | Low-loss ohmic RF-MEMS switches with interdigitated electrode topology[END_REF], [START_REF] Frlcke | Cantilever beam accelerometer based on surface micromachining technology[END_REF], [START_REF] Harness | Characteristic modes of electrostatic comb-drive X-Y microactuators[END_REF], [START_REF] Hoffmann | Electrostatic parallel-plate actuators with large deflections for use in optical moving-fibre switches[END_REF]. A wide variety of approaches has been proposed in literature to predict static and dynamic behaviour of microbeams [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF], [START_REF] Anac | Model validation methodology for design of micro systems[END_REF], [START_REF] Bona | Microsystems Mechanical Design[END_REF], [START_REF] Senturia | Microsystems design[END_REF]. Experimental validation is therefore aimed to verify their effectiveness in predicting pull-in condition and frequency response. Usually model sensitivity on the uncertainties of numerical values of design parameters and material properties is investigated. Very often is rather difficult to know precisely material properties and microspecimen dimensions.

FEM original approaches developed by the authors were proposed in [START_REF] Somà | Large deflections of microbeams under electrostatic loads[END_REF], [START_REF] Munteanu | Geometrical nonlinearities of electrostatically actuated microbeams[END_REF], [START_REF] Munteanu | A new FEM approach to the coupled-field analysis of electrostatic microactuators dynamics[END_REF], [START_REF] Brusa | Coupled-field FEM nonlinear dynamics analysis of continuous microsystems by non incremental approach[END_REF] and were already validated in [START_REF] Ballestra | Experimental characterization of the static behavior of microcantilevers electrostatically actuated[END_REF]. A preliminary experimental investigation was aimed to predict the static behaviour of planar microcantilevers, for inplane bending test. Present research is devoted to complete previous investigation activity focusing on out-of-plane bending microbeams.

Fully coupled electromechanical problem has electrical and mechanical coordinates, which are linked by the electromechanical coupling effect. Pull-in condition is responsible of a snap down of microbeam on the counterelectrode. In present case pull-in may be affected by some initial stress or strain present on the microsystem before the application of electric field. Moreover, it is well known that the problem is nonlinear because of the dependence of electromechanical force on displacement and voltage and, sometimes, of the so-called geometrical nonlinearity [START_REF] Munteanu | Geometrical nonlinearities of electrostatically actuated microbeams[END_REF], [START_REF] Munteanu | A new FEM approach to the coupled-field analysis of electrostatic microactuators dynamics[END_REF], [START_REF] Brusa | Coupled-field FEM nonlinear dynamics analysis of continuous microsystems by non incremental approach[END_REF], [START_REF] Ballestra | Experimental characterization of the static behavior of microcantilevers electrostatically actuated[END_REF]. To include these effects, models of microcantilever have to introduce the approximated analytical description of the initial curvature. It is usually sufficient to predict with enough accuracy the pull-in voltage and displacement, with an error of 2-3% maximum. For microbridges with double clamps axial stress is required to perform a coherent simulation of the actual system.

II. SPECIMEN CHARACTERIZATION

A. Fabrication process and measurement methods

Specimens used for this work were realized by ITC-IRST Research Center (Trento, Italy), by means of the so-called RF Switch (RFS) Surface Micromachining process. Gold is used for the suspended structures; material is deposited through electroplating by means of a chromium-gold PVD adhesion layer [START_REF] Margesin | Stress characterization of electroplated gold layers for low temperature surface micromachining[END_REF], [START_REF] Subhadeep | Stress and resistivity analysis of electrodeposited gold films for MEMS application[END_REF]. Profilometric measures and pull-in tests were performed by Fogale Zoomsurf 3D optical profiling system, based on non-contact optical interferometry [START_REF]Fogale Nanotech[END_REF]. The lateral resolution is ±0.3 µm, while the vertical resolution reaches ±0.5•10 -4 µm [START_REF] Ballestra | Experimental-numerical Comparison of Frequency Shift of Cantilever MEMS in Presence of Residual Stress Gradient[END_REF].Tables 1 and2 show the dimensions of microbeams used as specimens in testing, all measures are expressed in micron. In particular, t means beam thickness, L is length, w is width, g initial gap, y tip displacement for cantilever. Star symbol indicates measured value instead of nominal. Pull-in voltage (V PI ) was measured by gradually increasing voltage between suspended microbeam and ground, until the collapse of the structure on the ground counterelectrode.

Experimental pull-in was observed between the voltage range reported in tables I e II.

B. Residual stress

Presence of an initial curvature in microcantilever specimens is associated to an increase of pull-in condition. In double clamped microbeams pull-in depends on a pre-stress tensile loading caused by microfabrication process. Residual stress origin is supposed to be diffusion of chrome of adhesion layer among gold grains [START_REF] Subhadeep | Stress and resistivity analysis of electrodeposited gold films for MEMS application[END_REF] and difference of thermal expansion coefficient (CTE) between the gold structure and the photoresist layer underneath [START_REF] Stanec | Realization of low-stress Au cantilever beams[END_REF]. Residual stress grows every time there is a temperature difference during microfabrication process.

Stress is maximum at the interface gold/chromium/photoresist and decreases accross the microbeam thickness. In microcantilever, once removed the support substrate, traction at the interface is removed and turns into a deformation. Therefore initial curvature can be seen as an initial strain, which bends the microstructure out-of-plane [START_REF] Ballestra | Experimental-numerical Comparison of Frequency Shift of Cantilever MEMS in Presence of Residual Stress Gradient[END_REF]. In microbridges free bending is forbidden since microbeam is overconstrained by clamps. Therefore residual stress holds and strain is prevented .

III. MODELLING

According to the detailed descriptions of numerical models available in [START_REF] Rebeiz | RF-MEMS: Theory, design and technology[END_REF], [START_REF] Gaddi | Low-loss ohmic RF-MEMS switches with interdigitated electrode topology[END_REF], [START_REF] Frlcke | Cantilever beam accelerometer based on surface micromachining technology[END_REF], [START_REF] Harness | Characteristic modes of electrostatic comb-drive X-Y microactuators[END_REF], [START_REF] Hoffmann | Electrostatic parallel-plate actuators with large deflections for use in optical moving-fibre switches[END_REF], [START_REF] Anac | Model validation methodology for design of micro systems[END_REF], this paper will focus on the experimental validation of numerical results computed. Few modelling issues will be here resumed. Electromechanical force couples mechanical (displacements, rotations) and electrical (voltages, charges) degrees of freedom, and equilibrium equations in both the static and dynamic domains appear nonlinear. Moreover, typical dimensions of RF-MEMS introduce a second cause of nonlinearity, being referred to as "geometrical". In case of microcantilever a large displacement of the tip requires to resort to a nonlinear mechanical solution to find the actual equilibrium condition. This usually means to implement an iterative procedure which applies force step by step and locally linearizes the structural problem. In case of double clamped microbeam even in presence of small displacement the mechanical coupling occurring between the axial and flexural behaviour introduces either a hardening or softening effect on the structural stiffness. In both cases above mentioned, to solve the coupled problem a sequential solution is performed. Computation of the electric field distribution and of the related electromechanical force for a given equilibrium configuration of the deflected microbeam is separated from the mechanical analysis aimed to find the deformed shape of the beam under the electromechanical load. This sequence justifies implementing a computational loop. If the geometrical nonlinearity is active, a second iterative loop has to be implemented for each step of the electromechanical solution to find the actual equilibrium condition. Numerical methods are used to discretize the structure and the dielectric material. In present case both the dielectric and mechanical domains were meshed by FEM. A sequential approach based on Newton Raphson iteration method was implemented. The coupled electromechanical problem was solved by 2D and 3D models, by means of ANSYS code. In case of dynamic analysis Newmark Modified Algorithm was implemented and tested in connection with the non incremental approach for geometrically nonlinear structures. ANSYS code, MATLAB and FORTRAN implementations were performed and numerical results were compared to the experimental ones.

Alternately the coupled-field problem was solved by using coupled-field elements through a direct approach available in ANSYS as 1D transducer element TRANS126. It couples electro-mechanical domains and consists of a reduced-order model with structural translations and electric potential as degrees of freedom.

Initial deformed shape of microcantilver was analytically described by means of beam curvature κ and axial strain ε being written as function of flexural displacement v, axial displacement u, axial load N and bending moment M as follows:

;

2 2 EJ M dx v d = - = κ EA N dx du = = ε (1)
Generalized force vector F was then formulated as:

{ } { } [ ] [ ] [ ] ∑ ∫ ∑ ∫ = =               +               - = el el n i l T T T n i l T dx D B dx N M B F F 1 0 0 1 0 0 κ ε (2)
where initial values of force N 0 and Moment M 0 can be computed from the initial stress σ 0 ; ε 0T is the initial thermal strain and κ 0T the curvature in case of thermal contribution. To predict the actual pull-in of microbridge layout, axial effort had to be identified and included into the FEM models. This investigation was performed in ANSYS by applying either a distributed internal pressure at nodes or an initial strain as real constant. This action allowed estimating residual stress values by tuning the numerical pull-in tension on the experimental result.

IV. DISCUSSION

In more important and has to be evaluated by a 3D analysis to fit experimental result. In case of geometry 3 a prediction performed by ANSYS worse than by the non incremental approach was observed. Three dimensional effects are more effective in geometry 4, where force inputted in 2D models has to be tuned. Table 2 shows a large difference between pull-in voltage estimated without including residual stress and experimental measures. Microbridge 5 suffers the highest mismatch. Geometry 6 shows a good agreement and thus a pre-stress almost null. In case of geometry 7 pre-stress was present and detected. Tables demonstrate that upward initial curvature in geometrical nonlinear microcantilever is a very difficult configuration to be analysed. In fact, all 2D and 3D approaches implemented had to operate with a narrow gap and mesh morphing met some problems about pull-in, when the tip is close to the counterelectrode. Reduced order model based on TRANS126 gave better results, since it did mesh dielectric region. Double clamped microbeams strongly suffer pre-stress loading. A sensitivity analysis concerning initial tensile stress was performed. As Fig. 4 shows for geometry 5 different pre-stress conditions affect significantly pull-in voltage. It varies with axial stress from V PI = 29 (0 pre-stress) to 41 (10 MPa), 49 (20 MPa) and 57 (30 MPa). All models needed to be tuned by inputting a suitable value of tensile stress. In practice, all the microspecimens studied exhibited geometrical nonlinearity, thus requiring to resort to a nonlinear structural analysis based on the iterative solution of the mechanical problem within the sequential approach.

V. CONCLUSIONS An experimental investigation was performed to validate numerical approaches aimed to predict the behaviour of microbeams, nonlinearities due to electromechanical coupling and geometry were taken into account. The strong effect of residual stress on the pull-in voltage was detected and included in the models. The experimental results demonstrated a good agreement between the FEM approaches proposed, the methods allow a fairly fast and accurate prediction of the microbeams behaviour.

9- 11

 11 April 2008 ©EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5

Fig. 1 -

 1 Fig. 1-Geometry 4 series cantilever beams (L=200 um)

Fig. 2 .

 2 Fig. 2. Sequence of static equilibrium conditions measured by FogaleZoomsurf 3D during the pull-in tested performed on microcantilever.Geometry 2, sample 2.

Fig. 3 .

 3 Fig. 3. Sequence of static equilibrium conditions measured by FogaleZoomsurf 3D during the pull-in tested performed on microbridge. Geometry 5, sample 1.

Fig. 4 .

 4 Fig. 4. Analysis of the effect of the tensile load on geometry 5.

Table 1

 1

	calculated values of pull-in voltage of
	microcantilevers are compared to the measured ones. For
	geometry 1 all methods overestimated pull-in. In particular
	FEM 3D model was far from the true result. This was due to
	some problems of mesh morphing in presence of a very narrow
	gap. Geometry 2 shows a better agreement. Fringing effect is


T Veijola

Giorgio De Pasquale

A Somá

Comparison Between Damping Coefficients of Measured Perforated Micromechanical Test Structures and Compact Models

INTRODUCTION

Perforations are used in micromechanical squeeze-film dampers for several reasons. The main purpose is to reduce the damping and spring forces of oscillating structures due to the gas flow in small air gaps. Generally, the modeling problem is quite complicated, since the damping force acting on the moving structure depends on the 3D fluid flow in the perforations, in the air gap, and also around the structure. Compact models have been published in the literature, but their verification is generally questionable. Verification methods used are FEM solutions of the Navier-Stokes equations of the fluid volume and measurements [START_REF] Kim | Effect of Holes and Edges on the Squeeze Film Damping of Perforated Micromechanical Structures[END_REF].

In this paper, responses of four compact models are compared with measured responses of test structures. Six different perforated plates (figures 1 and 2) with different topologies have been measured at their first out-of-plane resonant frequencies, and the damping coefficients have been calculated from the quality factors (Q values) and effective masses. The measurement setup, the testing procedure and specimens characteristics are presented in [START_REF] Somà | Identification of Test Structures for Reduced Order Modeling of the Squeeze Film Damping in MEMS[END_REF]; here it is observed that dynamic parameters of the microsystem characterizing the fluidic and structural coupling can be extracted from the experimental frequency response function (FRF). The dynamic performance of microstructures are discussed based on the analytical solutions to perforated parallel-plate problems in [START_REF] Bao | Modified Reynolds' equation and analytical analysis of perforated structures[END_REF], [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF] and [START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF]. Since the perforation is uniform, the motion is almost translational, and also since the shape of the surface is rectangular, analytic damping models are applicable.

First, the oscillating flow is analyzed using several characteristic numbers, the applicable modeling method is then chosen, the damping coefficients are calculated and compared with the measured ones. Finally, the results are discussed.

II. TEST STRUCTURES AND MEASUREMENTS

Figures 1 and2 show the structures of the test specimen. The height of the plate h c = 15µm, the air gap height h = 1.6µm. Table I shows the other dimensions of the measured devices. In the table, q is the perforation ratio in percent, in this case q = MNs 0 2 /(LW), M and N are the number of holes in the length and width directions, respectively. The measurements are made using the interferometric microscope Fogale ZoomSurf 3D (figure 3), with 20x objective magnification factor, 0.1nm of vertical resolution and 0.6µm of lateral resolution.

The Frequency Shift technique is used, consisting of the excitation of the structure by an alternate voltage, the frequency of which is progressively increased by discrete steps. For each level of actuation frequency, the corresponding amplitude of vibration is stored and the experimental FRF is plotted for the detection of the resonance peak (figure 4). The first detection is made across a wide frequency range (0-500kHz) in order to roughly locate the resonance; five successive identical detections are then performed across a more precise and narrow range. These are statistically treated to extract the values of resonance reported in Table II. The measurement technique described uses a red monochromatic light source for the interferometric fringes detection. The vibration amplitude is detected optically in correspondence of a selectable region (detection window) of the specimen, located at the center of the suspended plate. The output value of the oscillation amplitude is averaged between the values captured by each pixel of the CCD camera inside the active window. The quality factor is extracted from the experimental curve, that was previously interpolated by a 6-order polynomial; the damping coefficient is finally calculated from the quality factor, resonant frequency and the effective mass by means of the method of the half power bandwidth. These are shown in Table II. The effective mass is calculated from FEM eigenmode analysis. The mass ratio α is the ratio between the modal mass and total mass. It is ensured that the amplitude of the oscillation is small compared with the gap height, and the static bias voltage caused by the excitation signal and DC bias voltage does not deflect the plate changing the air gap height.

III. MODELING OF THE DAMPING COEFFICIENT

A. Analysis of the fluid flow

The actual mass is supported with thin beams in such a way that the movement of the mass is approximately translational. This justifies starting with an analysis where the velocity of the plate surface is constant.

Flow patterns

In perforated dampers in perpendicular motion, two different flow patterns can be distinguished. The first is the "closed borders" pattern, where the fluid flows only through the holes. The second pattern, the "closed holes" pattern considers only the flow from the damper borders. In practical dampers both patterns exist simultaneously. The perforation ratios q (the area of the surface without holes divided by the area of the holes) are here considerably high, ranging from 24% -59%. Also, the holes are relatively wide compared to the air gap height. It is then evident that the first "closed borders" flow pattern is strong here. Here compact models that consider both flow patterns are selected; the contribution the flow patterns in the measured cases will be discussed later in this paper.

Rarefied gas

Air at standard atmospheric conditions is used in the measurements. The pressure P A = 101kPa, the density ρ = 1.155 kg/m 3 , the viscosity coefficient µ = 18.5⋅10 -6 Ns/m 2 , and the mean free path λ = 65nm. The air gap height h = 1.6µm, which makes the Knudsen number of the air gap flow K ch = λ/h to be 0.04, and the smallest hole diameter is 5µm making the Knudsen number of the perforation K tb = λ/s 0 to be 0.013. The estimated contributions to the damping coefficient are 1/(1 + 6K ch ) and 1/(1 + 7.567K tb ), that is -24% and -9.8%, respectively. The assumption of a slightly rarefied gas is justified, and the slip velocity model is sufficient.

Compressibility

Next, the contribution of compressibility is analyzed. This can be made by studying the squeeze number. For a rigid surface, the squeeze number is W is here the smallest, "dominating" characteristic dimension. For example, for case A without holes ω σ

6 10 8 . 3 - ⋅ =
(2) At 200kHz σ = 4.8. This means that without perforation, compressibility should be considered (when σ ∼ 20, the viscous and spring forces are equal). When the surface is perforated, the situation changes completely. The characteristic dimension can now be estimated to be the space between the holes. In case A this is 5.2µm. The squeeze number now becomes σ = 0.03 at 200kHz. This is an overestimate of the real situation, since accounting for the damping in the holes will make the effective squeeze number considerably smaller.

According to the squeeze number analysis, the spring forces are negligible compared with the damping forces, and noncompressible gas can be assumed without loss of any accuracy. The damping coefficient can be considered constant up to several MHz. Also, it is expected that the spring force due to the gas is much smaller than the force due to the effective spring of the mechanical structure. The frequency shift due to gas compressibility is expected to be very small. The spring force in the system is only due the effective spring of the massspring system. The spring coefficient can also be considered constant at least up to several MHz.

Gas Inertia

One can suspect that the inertia of the gas may contribute to the damping coefficient. The place where the inertia is the largest is the "widest" flow channel, that is the perforation holes.

The contribution of inertia is characterized by the Reynolds number R e , specified for a circular channel as [START_REF] Morris | Oscillatory Flow in Microchannels[END_REF]:

µ ω ρ 2 e r R = (3)
The "worst" case, where the inertial is the largest, is case D where the hole diameter is the largest. If the square channel is approximated with a circular channel having a radius of r = 4µm, that's half of the hole side, gives R e = 0.998•10 -6 ω, and at 200 kHz R e = 1.255. The real and imaginary parts of the impedance are equal when R e ~ 6. The additional imaginary part does not directly influence the damping coefficient, but the change is due to the frequency dependent real part. At R e ~ 6, the change in the real part in only 3.2% [START_REF] Morris | Oscillatory Flow in Microchannels[END_REF].

This Reynolds number study shows that the inertia needs not to be considered even in the accurate analysis.

B. Compact models

A model for the noncompressible perforation cell is sufficient, as indicated in the study above. Four models that consider both "closed holes" and "closed borders" flow patterns are selected to be compared. The first one, M1, is a model by Bao [START_REF] Bao | Modified Reynolds' equation and analytical analysis of perforated structures[END_REF] for a rectangular damper that is much longer than wide. In this model the air gap regime flow resistance, the flow resistance of a circular perforation, and its constant elongation are included. Continuum flow conditions are assumed. The 2 nd model M2 has been also presented in [START_REF] Bao | Modified Reynolds' equation and analytical analysis of perforated structures[END_REF], but now an arbitrary rectangular surface is assumed. Next, the model M3 in [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF] is used. The air gap flow resistance model, the circular perforation flow channel model, and four different elongations of the flow channels, that vary depending on the ratios of the cell dimensions, are included in the model. Slip velocity conditions are used for the air gap and the perforations. The 4th model, M4, is made especially for square holes [START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF]. It includes similar components as the previous model and accounts also for the rarefied gas in the slip flow regime.

Model M4 for square perforations accepts directly the dimensions given in Table I (size of the perforation cell s x = s 0 + s 1 ). To apply the other models, an effective radius of the circular perforation r 0 and the cell r x need to be specified first. Matching the areas of the actual cell s x 2 , and the equivalent circular cell πr 0 2 gives

π x x s r = (4)
The radius r 0 is determined by requiring the acoustic impedances of square and rectangular channels to match. For relatively small Knudsen numbers this leads to approximately [START_REF] Bao | Modified Reynolds' equation and analytical analysis of perforated structures[END_REF], [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF] The Appendix shows all equations needed in computing the damping coefficients using models M1…M4.

IV. SIMULATION RESULTS AND DISCUSSION

The results of the comparison are shown in Table III. ∆ i is the relative error of the simulated damping coefficient c s of model Mi compared to the measured damping coefficients. The results of models M3 and M4 are quite close to each other, the largest error between them is only 2.8%-points, showing that the effective radius approach is sufficient. Continuum conditions were assumed for M1 and M2, giving approximately 10% larger values than with slip velocity conditions. If these models are corrected to account for the III.

The drag on the sidewalls is expected to increase the damping coefficient since the structure is relatively high. The models do not account for this drag force. The moving supporting beams will have also an additional contribution to the damping. The length of the supporting beams in all cases is about L b = 122µm, and their widths are about W b = 4µm. A rough approximation for the contribution of the supporting beams is

) 6 1 ( ) 3 . 1 ( 3 4 ch 3 3 b b b K h h W L c + + = µ (6)
This gives a damping coefficient of c b = 0.16⋅10 -6 Ns/m. This approximation shows that the damping due to the beams is very small.

The responses of models M1 and M2 are quite close. One could expect that the error of M1 would be larger in case E and especially in case F, since the length to width ratio is quite small in these cases. The explanation for this can be found by studying the contribution of the different flow patterns. This can be easily done using the perforation cell model that assumes the "closed borders" flow pattern, where the pressure distribution is independent of the shape of the damper. The flow resistances R P for a perforation cell in [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF] and [START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF] are derived using this assumption. The damping coefficient becomes simply c P = NMR P (MN is the number of holes). Table IV shows the errors of the damping coefficients c P compared with the measured values using the models for circular cells, M5, and rectangular cells, M6, as presented in [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF] and [START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF], respectively. The results in Table IV show that the "closed borders" flow pattern is the dominant one. The contribution of the "closed holes" flow is only 6% -16% of the damping coefficient. This explains why models M1 and M2 differ only slightly in this case: the contribution of the shape-dependent damping is quite small. To study further the sources of damping, Table V shows the contributions of the flow resistance components in M5 [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF]. The flow resistance of the perforations R S is the most significant source of damping; its contribution is approximately 65%. The 2 nd important contribution comes from the intermediate region resistances R IS , R IB , and R IC : 15 -20%. Next important is the elongation at the perforation outlet R E , about 8%.

V. CONCLUSIONS

Measured damping coefficients have been compared to those obtained with four different compact models for perforated dampers. After analyzing the oscillating flow with several characteristic numbers, sufficient models were selected. Only translational motion was assumed. The results of all models were quite close to each other, a systematic underestimate of the damping coefficient was about -20%. The reasons for this were discussed and the contribution of various flow components were presented. For a more accurate analysis, the realistic modes of the plates should be considered. It is expected that the "closed holes" flow pattern will become relatively stronger in this case.

The comparison showed also how a model for circular perforations can be used to model square holes.

ISBN: 978-2-35500-006-5 APPENDIX This Appendix contains equations for four compact models M1…M4. The dimensions and symbols in Fig. 1 are used: the length and width of the perforated plate are L and W. The side lengths of the square holes and the square perforation cells are s 0 and s X = s 0 + s 1 , respectively.

A. Model M1 equations

The equations for a narrow hole plate (L>>W) are given in [START_REF] Bao | Modified Reynolds' equation and analytical analysis of perforated structures[END_REF]. Note, in the following equations a = W/2 and b = L/2. The equivalent radii for the circular cell and hole are given in Eqs. ( 4) and [START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF]. The damping coefficient c is

            -         + = l a a l h h K r r h aL c C C tanh 1 16 ) ( 3

B. Model M2 equations

The equations for an arbitrary shaped rectangular plate are also included in [START_REF] Bao | Modified Reynolds' equation and analytical analysis of perforated structures[END_REF]. Also, in the following equations a = W/2 and b = L/2. The equivalent radii for the circular cell and hole are given in Eqs. ( 4) and [START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF]. The damping coefficient c is ( ) ( )

C. Model M3 equations

A model for a circular perforation cell is derived in [START_REF] Veijola | Analytic Damping Model for an MEM Perforation Cell[END_REF], and the damping coefficient of a rectangular perforated plate is given in the paper. Note, in the following equations a = W and b = L. The equivalent radii for the circular cell and hole are given in Eqs. ( 4) and ( 5). The damping coefficient c is

m MNR R Q h n m b n a m b a G n m n m =         + =
The flow resistance of a single perforation cell is

( ) ( )         ∆ + = + ∆ = ∆ = ∆ - =         - + - = + + + + + = 0 E tb C E C C 0 IC B 0 IB S 2 0 2 2 0 2 X IS 4 X 4 0 2 X 2 0 0 X 3 ch 4 X S E C IC 4 0 4 X IB IS S P 8 8 8 6 8 2 8 3 ln 2 1 12 r Q h R R r R r R h r r r R r r r r r r h Q r R R R R r r R R R R πµ πµ πµ πµ πµ
where the elongations are

      + +         - = ∆ + + - = ∆ h h h r f K K r r K r r r r C 0 B ch tb 2 X 2 0 B ch 2 X 2 0 X 0 S , 1 732 . 0 1 812 . 0 1 33 . 1 5 . 2 1 86 . 0 32 . 0 56 . 0 ( ) ( )               - + × + ⋅ = ∆         - - + = ∆ h r f r r r r K r r r r K 0 E 4 X 4 0 2 X 2 0 tb E 2 X

Fig. 1 .

 1 Fig. 1. Geometrical shape and dimensions of the vibrating structures.

9- 11

 11 April 2008 ©EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5

Fig. 3 .Fig. 4 .

 34 Fig. 3. Interferometric microscope Fogale ZoomSurf 3D (a) and the 20x Nikon objective (b).

9- 11

 11 April 2008 ©EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5

9- 11

 11 April 2008 ©EDA Publishing/DTIP 2008 ISBN: 978-2-35500-006-5 rarefied gas, the errors become approximately 10% worse, than those shown in Table

  is the same as used in M1 equations.

  The flow rate coefficients and Knudsen numbers for the air gap and the holes are 9A model for a rectangular perforation cell has been given in[START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF]. Note, in the following equations a = W and b = L. The damping coefficient c is given by (C1), where R P for a rectangular hole isThe equation for ∆ E includes a misprint in[START_REF] Veijola | Analytic Damping Model for a Square Perforation Cell[END_REF]. The corrected equation is shown above.The flow rate coefficients and Knudsen numbers for the square