A numerical study of squeeze-film damping in MEMS-based structures including rarefaction effects

In a variety of MEMS applications, the thin film of fluid responsible of squeeze-film damping results to be rarefied and, thus, not suitable to be modeled though the classical Navier-Stokes equation. The simplest way to consider fluid rarefaction is the introduction of a slight modification into its ordinary formulation, by substituting the standard fluid viscosity with an effective viscosity term. In the present paper, some squeeze-film damping problems of both parallel and torsion plates at decreasing pressure are studied by numerical solving a full 3D Navier-Stokes equation, where the effective viscosity is computed according to proper expressions already included in the literature. Furthermore, the same expressions for the effective viscosity are implemented within known analytical models, still derived from the Navier-Stokes equation. In all the considered cases, the numerical results are shown to be very promising, providing comparable or even better agreement with the experimental data than the corresponding analytical results, even at low air pressure. Thus, unlike what is usually agreed in the literature, the effective viscosity approach can be efficiently applied at low pressure regimes, especially when this is combined with a finite element analysis (FEA). SOMMARIO. In molte applicazioni MEMS, il sottile strato di fluido responsabile del fenomeno dello squeezefilm damping risulta essere rarefatto e, quindi, non modellabile mediante l’equazione classica di Navier-Stokes. Il modo più semplice per tener conto della rarefattezza del fluido consiste nell’introduzione di una piccola modifica all’interno della sua formulazione ordinaria, ovvero nella sostituzione della viscosità del fluido con una viscosità effettiva. Nel presente lavoro, vengono presi in considerazione alcuni problemi di squeeze-film damping, riguardanti casi in cui le piastre coinvolte sono dotate di moto normale traslatorio e casi in cui esse sono dotate, invece, di moto torsionale. Tali problemi vengono risolti mediante una modellazione numerica 3D dell’equazione di Navier-Stokes, in cui la viscosità effettiva viene calcolata mediante apposite espressioni già note in letterature. Inoltre, queste stesse espressioni sono utilizzate all’interno di modelli analitici già conosciuti, derivati anch’essi dall’equazione di Navier-Stokes. In tutti i casi qui considerati, i risultati numerici sono molto promettenti, in quanto essi sono caratterizzati da uno scostamento dai dati sperimentali uguale o inferiore a quello fornito dai corrispettivi risultati analitici, anche a basse pressioni. Quindi, nonostante quanto sia usualmente asserito in letteratura, il metodo della viscosità effettiva può essere efficientemente adottato anche a bassi regimi di pressione, specie se in combinazione con analisi numeriche agli elementi finiti (FEA).


INTRODUCTION
queeze-film damping is the main source of energy dissipation affecting MEMS devices, like microaccelerometers, microgyroscopes and micromirrors.It is related to the presence of a thin film of fluid, confined between a solid moving surface and a substrate.The reciprocal movement of the two surfaces causes the fluid to be sucked into/pulled out of its tight channel.This generates a pressure field inside the fluid, which is then responsible of a resistive force to hinder the walls' movement.Squeeze-film damping becomes significant when the thickness of the fluid layer is equal or smaller than one third of the width of the confining surfaces [1].Such condition is accomplished in many MEMS applications, and it is the reason why studies about squeeze-film damping have increased since the 1990s.Generally, the fluid flow through its tight channel is modeled through the classical Navier-Stokes equation, which simplifies in the Reynolds equation under some assumptions (e.g., inertial effects, thermal gradients, and fluid out-of-plane movements are negligible).However, the main hypothesis for the Navier-Stokes/Reynolds equation to be valid is the continuity of the fluid body.Such approximation applies at high pressure (e.g., atmospheric pressure), but fails at low pressure.Since MEMS-based devices may work in vacuum, it is important to have accurate modeling of squeeze-film damping even in this condition.The parameter conventionally used to evaluate fluid rarefaction is the Knudsen number, K n .This is defined as the ratio of the mean free path of the fluid molecules (λ) to the characteristic dimension of the fluid channel (l): where λ is [2]: where R is the gas constant, T is the temperature, d is the diameter of the fluid molecules, NA is the Avogadro number, and p is the ambient pressure.For low values of K n (K n <0.01: continuum regime), the Navier-Stokes equation is valid, whereas for medium (0.01<K n <1: slip regime/transition regime) and high (Kn>1: free molecular regime) values of Kn, continuum approximation does not apply and other models than the Navier-Stokes equation should be considered [1].
In spite of the intense efforts already done by the scientific community on modelling fluid rarefaction, there is still need of further improvement and understanding.
In the present paper, a numerical approach is adopted to study some squeeze-film damping problems, involving either normal or torsional movement of the moving plate at different pressures, ranging from the atmospheric pressure to almost vacuum.The numerical results are compared to the experimental data, already available in the literature, and to the results obtained by implementation of known analytical models, derived from the Navier-Stokes equation.

MODELING OF FLUID RAREFACTION
n the literature two main approaches can be found for modeling of fluid rarefaction.The first one is based on the introduction of a slight modification in the Navier-Stokes equation.In particular, the standard fluid viscosity (μ), which compares into its classical formulation, is substituted with a scaled quantity, known as effective viscosity (μ eff ), defined as: Where QPR is the flow rate coefficient, computed as a function of the Knudsen number.During the years, many expressions for Q PR , and consequently for μ eff , have been proposed.However, in the present paper, the attention is focused on those, which have been proved to work better [3][4].Thus, three expressions are considered herein, which were proposed in [5], [6], and [4], respectively.S I According to Veijola et al. (1995), the effective viscosity (μeV) can be computed as [5]: The effective viscosity (μ eL ) proposed by Li (1999) is instead [6]: where a=0.01807, b=1.35355, and c=-1.17468, and D is the inverse Knudsen number, defined as: Pandey and Pratap (2008) further improved Li's model, in order to achieve better agreement with experiments.Thus, they introduced a scaling factor for D equal to 1.4, and proposed their own expression (μ eP ) as [4]: where a, b, c, and D have the same meaning as before.
Since its simplicity, the effective viscosity approach could be very powerful.However, none of the above expressions have been proved to work well in all conditions [7].
The second approach for modeling fluid rarefaction is based on the molecular dynamics (MD).As opposite to the first approach, in this case attention is paid to collision of the fluid molecules with the walls' surface, whereas interactions between fluid molecules are neglected.In 1966, Christian was the first to adopt the MD for computation of the quality factor (Qc) related to squeeze-film damping of a plate moving normal to the substrate.The expression he proposed was [8]: where ρ is the density of the solid surface, b the beam thickness, fr the resonance frequency, R the gas constant, T the temperature, M the molar mass of the fluid molecules, and P the ambient pressure.Such model was validated through some experiments on miniaturized beams, performed in [9].Because of the resulting poor agreement with the experimental data, this model was further improved by Bao et al. [10], who proposed their own expression (QB): where L is the peripheral length of the plate confining the fluid, and d0 is the initial thickness of the fluid layer.With respect to Christian's model, which is derived from momentum transfer between fluid molecules and the solid moving surface, the model proposed by Bao et al. is based on energy transfer, and it allows for consideration of the size of the surface confining the fluid and the effects due to the presence of other fluid in the surroundings.Such model was further improved by Hutcherson and Ye [11], who performed MD simulations in order to relax some constraints of the previous model, like constant particle velocity and constant beam position.However, all the aforementioned MD models were suitable for describing squeeze-film damping in case of rigid structures.Only recently, Li and Fang [12] improved Hutcherson and Ye's model for description of even torsion and flexible beams.
Since the models based on MD do not consider interactions between fluid particles, they do not take into account the viscous character of the fluid.Thus, they are suitable to describe squeeze-film damping when the fluid is highly rarefied, e.g. the free molecular regime.Furthermore, the simplified equations reported herein have been proved to not provide sufficient agreement with experiments [3,9] and the more recent MD based models require the development of proper more or less laborious numerical codes.The first approach, being instead based on the Navier-Stokes equation, considers only the viscous contribution to squeeze-film damping.Thus, it should be effective only in the continuum and transition regimes.However, because of its simplicity it was also adopted in the free molecular regime, providing reasonable errors [3][4], too.The way such approach is usually implemented for computation of squeeze-film damping in terms of damping coefficient (or quality factor) requires the substitution of the standard viscosity with an effective viscosity term in the Navier-Stokes equation or in compact formulae derived from that.However, such formulae are available only for regular and simple geometries, where the plate confining the fluid has either parallel or torsional movement with respect to the substrate [1].Besides, it was proved by the authors of the present paper that at high pressure solution of the Navier-Stokes equation by numerical analysis can be more effective than simplified formulae [13].Thus, in the following, the effective viscosity approach will be adopted, coupled with numerical solution of the Navier-Stokes equation.

ANALYTICAL MODELING OF SQUEEZE-FILM DAMPING IN RAREFIED REGIME
he literature provides compact formulae, derived from the solution of the Navier-Stokes equation, to describe squeeze-film damping in rigid MEMS structures, moving either normal or torsional with respect to the substrate.When the thin film of fluid is confined between a substrate and a rectangular plate, which moves normal to the substrate, the damping coefficient can be determined as [1]: where L and w are the plate length and width, respectively, h is the thickness of the fluid film, β is a correction factor, depending on the w/L ratio, and μ is the fluid viscosity, which can be substituted with the effective viscosity, computed according to one of expressions ( 4), ( 5), and ( 6).
There is an alternative semianalytical formula, which is valid in case of a rigid plate, moving normal to the substrate.This is [14]: where ω is the frequency of the plate movement, and QPR, Gmn, and Cmn are defined as:     where h and μ are the same as before; / q j   , ρ is the fluid density, λ is the mean free path, p is the ambient pressure, a and b are the plate sides, and n γ is a coefficient depending on heat conduction and temperature boundary conditions.As opposite to compact formula (9), here it is not possible to isolate an effective viscosity term, which could be changed according to other expressions.In fact, rarefaction terms are embedded within the whole expression.When the plate confining the fluid is provided with torsional movement, the compact formula to be used to determine the damping coefficient is [15]: where h, μ, L, and w are the same as before, and η=w/L.Similarly to Eq. ( 9), in order to take into account fluid rarefaction, the fluid viscosity μ can be substituted with the effective viscosity, computed according to (4), ( 5), or (6).

NUMERICAL MODELING OF SQUEEZE-FILM DAMPING IN RAREFIED REGIME
ince the versatility and computational power of modern computers, a numerical approach for solving squeeze-film damping problems can be a valid alternative to analytical formulae, especially when complex geometries have to be studied.In the present paper, numerical analyses of squeeze-film damping problems were performed by the use of a commercial finite element software, Comsol Multiphysics, at different pressure regimes, ranging from the atmospheric T S pressure to almost vacuum.The main advantage of such software is its capability to perform multi physics analysis.Thus, it is particularly suitable to deal with problems defined in different physical domains, like these, which contemporarily involve structural mechanics and fluid dynamics.The effective viscosity approach was adopted to take into account fluid rarefaction.Thus, a full 3D incompressible Navier-Stokes equation was solved, which included both the air between the moving plate and the substrate, and the air in their surroundings.Rarefaction of the air between the moving plate and the substrate, was modeled by computing the effective viscosity through one of expressions ( 4), ( 5), or (6).To model rarefaction of the air in the surroundings of the moving plate, the effective viscosity approach cannot be adopted, since it consists of scaling the standard fluid viscosity by a factor depending on the Knudsen number.Nevertheless, the Knudsen number depends on the characteristic length of the channel, where the fluid molecules flow, and such channel is very wide for the air in the working volume, causing the correction factor to be almost one.Thus, the viscosity was lowered according to another procedure, based on the following considerations.The mean free path of the fluid molecules is inversely proportional to the ambient pressure at constant temperature (Eq.( 2)).It is instead proportional to the fluid viscosity, according to the following equation [16]: where R is the individual gas constant, which is equal to 286.9 JK -1 kg -1 for air.In the present numerical analyses, the mean free path could not be increased to effectively simulate fluid rarefaction.Thus, the viscosity was reduced according to expressions ( 2) and ( 15) at low pressures (when the Knudsen number inside the fluid channel is larger than 1, corresponding to the free molecular regime), while keeping the mean free path constant.The software automatically generates a 3D mesh, consisting of tetrahedral elements modeling the moving plate, the fluid under the plate, and the fluid in its surroundings.In particular, the volume of fluid to consider in the analysis was found to not affect the results, if it extends to a region sufficiently far from the plate edge, in order for the fluid flux to develop completely.Fig. 1 shows a typical mesh generation and the pressure field under the moving plate.The number of mesh elements ranged from about 5000 to 50000, depending on the particular case, and was decided after a convergence study.In all the analyses reported herein, geometric symmetry was taken into account, in order to reduce the computational work, which was performed by a workstation with the following technical features: RAM 16 GB, Intel(R) CORE(TM) i7 CPU 860 @ 2.80 GHz.In these conditions, the time required for one simulation ranged from 15 to 60 minutes.Two kinds of squeeze-film damping problems were considered, each involving either a plate moving normal or torsional with respect to the substrate (Fig. 2).
For the problems where the plate moves normally, the damping coefficient was evaluated as: where p(x,y) is the pressure field over the moving plate and v z its velocity.In the second kind of problems, where the plate rotates around the x-axis (Fig. 2), the damping coefficient was evaluated as [17]: where ω=2πf, being f the frequency of the plate movement.

PLATES
o evaluate the effectiveness of the numerical analysis, the squeeze-film damping problems reported by Sumali [3] are considered herein.He investigated a gold non-perforated plate (Fig. 3), moving normal to the substrate, and determined the corresponding damping coefficient at different pressures, ranging from the atmospheric value to few Pa.For such geometry, three sets of numerical analyses were performed [18], each implementing one of expressions ( 4), ( 5), and ( 6) for the effective viscosity of the fluid inside the channel.Fig. 4a shows a log-log plot reporting a comparison between the numerical results (c n ) and the experimental data (c e ) at varying pressures.Because of its complex profile, analytical expressions ( 9) and ( 10) are not suitable for the geometry in Fig. 3. Thus, these were applied to study an equivalent rectangular plate, with the same area as the original surface, whose size is given in [3].Furthermore, to take into account fluid rarefaction in Eq. ( 9), the standard fluid viscosity was substituted with the effective viscosity, obtaining three sets of data (one for each expression).Fig. 4b shows a log-log plot reporting a comparison between the analytical results (ca) and the experimental data (ce), as functions of pressure.Both the analytical and the numerical approach work well from high to medium pressures (around 103 Pa).At low pressure instead, the numerical results show much better agreement with the experimental data (average difference of 23%) than the results obtained with Eq. ( 9) (average difference of 58%).This can be explained since Eq. ( 9) does not take into account border effects, which are instead considered in the numerical simulations.In order to overcome such limitation of the analytical model, a correction factor should be introduced.The only elongation model, providing a correction factor, available in the literature, was verified for Kn smaller than 0.13.However, most of the experimental points reported in [3] do not respect such condition.Thus, the corresponding correction factor, being not valid, was not included in the graph of Fig. 4b (and in the following graphs), even if that would improve the results.Eq. ( 10) shows better agreement with experiments than numerical analysis at low pressure, but it does not provide an expression for the effective viscosity to be easily implemented in numerical simulations.Thus, a direct comparison is not possible.

MICROMIRRORS
n this section, squeeze-film damping affecting three square torsion micromirrors was considered at varying pressures, ranging from the atmospheric pressure to almost vacuum.Such micromirrors were experimentally investigated by Minikes et al [17] and Pandey and Pratap [4], and their geometry and resonance frequency are reported in Tab.

Minikes et al. (2005)
provided their experimental results in terms of quality factor (Q).This is related to the damping coefficient (c) as: Being Ix the plate inertia moment around the rotation axis x.Figs.5a and 5b shows how the quality factor varies with pressure in both micromirrors.At very low pressures (below 20Pa), the quality factor becomes almost constant.This occurs since in such conditions squeeze-film damping is no longer the main damping component, but structural damping prevails.Thus, in order to correctly evaluate squeeze-film damping at such low pressures, the structural component has to be isolated.This can be done according to the following procedure [3].
A restricted number of experimental points, each corresponding to very low pressures, is considered and interpolated by linear data fitting.The intersection point of the interpolation curve with the y-axis (e.g.quality factor at 0 Pa) corresponds to the quality factor associated to the structural damping (Fig. 6).The value of the structural quality factor are 15886 and 5366 for the two cases, respectively.For the cases considered herein, the quality factor was determined by both the analytical and numerical approach.In particular, in Eq. ( 14) the standard viscosity was substituted with the effective viscosity, computed by expressions (4), ( 5), and (6), obtaining three sets of analytical data, respectively.Similarly, in the numerical analysis three sets of simulations were performed, each implementing one expression for the effective viscosity.In order to take into account the structural damping, a total quality factor (Q tot ) was computed as [4]: where Qsq is the quality factor computed according to either the numerical analysis or the analytical formula, and Qst is the structural quality factor.Fig. 7 reports a log-log plot showing a comparison between the numerical results [20] and the experimental data for the two torsion mirrors, investigated by Minikes et al. [17] at varying pressure.Similarly, Fig. 8 reports a comparison between the analytical results and the experimental data.From such comparison, it is possible to notice the good agreement with experiments provided by both the numerical analysis and the analytical formula at all pressures.In this case, the difference between the numerical and the analytical results is less significant than in the previous case.In fact, the average difference between the numerical data and the experimental ones is 24%; while the analytical model provides an average difference of 27%.The reason why both the numerical and the analytical modeling offer similar results could be related to the border effects, which in these cases do not play a major role on damping.
The final problem considered herein is the torsion micromirror reported by Pandey and Pratap [4].In this case, it is not necessary to compute the structural damping, since the authors already provided corrected experimental data, which do not take it into account.Similarly to the previous case studies, both the analytical and the numerical approach were applied to determine squeezefilm damping in terms of quality factor.As before, three sets of numerical results and analytical results were obtained and compared to the experimental data in Fig. 9a and 9b, respectively.In this case, the difference between the analytical and the numerical results is much more significant than before (the average difference between the numerical and the corresponding experimental data is 18%, while it can be even an order of magnitude for the analytical results).This is related to the geometry of the present problem, where the ratio of the fluid gap thickness to the plate width is significantly larger, and accordingly the magnitude of the border effects, which are neglected when applying formula (14).However, also in this case the effectiveness of the numerical analysis emerges, providing results in very good agreement with the experiments.

CONCLUSIONS
n the present paper, four different squeeze-film damping problems were considered, at varying pressures, ranging from the atmospheric value to almost vacuum.To analyze such cases, involving both normal and torsion movement of the plate confining a thin film of air, both an analytical and a numerical approach were adopted, both of them based on the Navier-Stokes equation.In order to model fluid rarefaction, the effective viscosity approach was followed, which consists of substituting the standard fluid viscosity, contained in the equation, with a scaled term, known as effective viscosity.The literature provides many expressions for computing the effective viscosity, and the three expressions, which were proved to work better, were considered within the paper.Then, each of the four considered problems was solved by both the numerical and analytical methods, implementing each time one of those expressions, obtaining four sets of numerical results and four sets of analytical data.In all the considered cases, the numerical results were very promising even at very low pressure, with values of the squeeze-film damping coefficient/quality that were comparable (and even closer than the analytical results) to the experimental data.Thus, unlike what is usually agreed in the literature, the effective viscosity model, especially when combined with FEM analysis, results to be effective in a wide range of pressures, including very low values.In addition, thanks to the computational power of modern computers and the versatility of FEM analysis, the procedure described herein can be implemented in a variety of applications, extending to complex and more realistic structures, as opposite to analytical approaches.In this work, the effective viscosity was computed according to already known expressions.However, a future work can be focused on derivation of a new expression to be implemented in the numerical analysis, in order to get results even closer to the experimental data.

Figure 1 :
Figure 1: Mesh of tetrahedral elements (a) and pressure field under one fourth of the moving plate (b).

Figure 2 :
Figure 2: Reference system for the moving plate.

Figure 5 :
Figure 5: Values of the damping coefficient versus pressure of the torsion micromirrors 1 (a) and 2 (b) reported by Minikes et al. [17].

Figure 6 :
Figure 6: Procedure to determine the quality factor associated to the structural damping for micromirror 1 (a) and 2 (b) reported by Minikes et al [17].

Figure 7 :
Figure 7: Comparison of the numerical results with the experimental data for micromirror 1(a) and 2 (b) reported by Minikes et al [17].

Figure 8 :
Comparison of the analytical results with the experimental data for micromirror 1 (a) and 2 (b) reported by Minikes et al[17].

Figure 9 :
Figure 9: Comparison of the numerical (a) and the analytical results (b) with the experimental data reported by Pandey and Pratap [4].