A Monolithic, FEM-Based Approach for the Coupled Squeeze Film Problem of an Oscillating Elastic Micro-Plate Using 3D 27-Node Elements

In this study we describe an FEM-based methodology to 
solve the coupled fluid-structure problem due to squeeze film effects present 
in vibratory MEMS devices, such as resonators, gyroscopes, and acoustic 
transducers. The aforementioned devices often consist of a plate-like 
structure that vibrates normal to a fixed substrate, and is generally not perfectly 
vacuum packed. This results in a thin film of air being sandwiched between the 
moving plate and the fixed substrate, which behaves like a squeeze film 
offering both stiffness and damping. Typically, such structures are actuated 
electro-statically, necessitating the thin air gap for improving the efficiency 
of actuation and the sensitivity of detection. To accurately model these 
devices the squeeze film effect must be incorporated. Extensive literature is 
present on mod- eling squeeze film effects for rigid motion for both perforated as 
well as non-perforated plates. Studies which model the plate elasticity often 
use approximate mode shapes as input to the 2D Reynolds Equation. Recent works 
which try to solve the coupled fluid elasticity problem, report iterative 
FEM-based solution strategies for the 2D Reynolds Equation coupled with the 3D 
elasticity Equation. In this work we present a FEM-based single step solution 
for the coupled problem at hand, using only one type of element (27 node 3D 
brick). The structure is modeled with 27 node brick elements of which the 
lowest layer of nodes is also treated as the fluid domain (2D) and the integrals 
over fluid domain are evaluated for these nodes only. We also apply an 
electrostatic loading to our model by considering an equivalent electro-static 
pressure load on the top surface of the structure. Thus we solve the coupled 
2D-fluid-3D-structure problem in a single step, using only one element type. 
The FEM results show good agreement with both existing analytical solutions and 
published experimental data.

The wide sca MEMS senso to increasing due to the thin cally such de brating norma trapped in-bet sions of the height of the spring and a squeeze film nant dissipatio operating in * Corresponding a plied Mathemati ne November 20 g/10.4236/jamp.

Monolith ze Film
Anish 1  the Reynolds Equation [3].With rigid plate assumption the Reynolds Equation can be decoupled from the elasticity Equation and, further, on linearization can be solved to obtain analytical expressions for stiffness and damping.Blech [3] studied the effect of squeeze film induced stiffness and damping for rigid plates with trivial pressure boundary conditions.Darling et al. [4] presented analytical solutions to the linearized Reynolds Equation for various venting conditions, using a Greens function approach.For flexible structures one has to account for the variable air gap, and the elasticity Equation has to be coupled with the Reynolds Equation for accurate modeling.Hung et al. [5] presented a reduced order macro model based on basis functions generated from finite difference simulations.They applied this technique to model a pressure sensor as a clamped-clamped microbeam to study the pull-in dynamics of the system using a 1D Euler beam Equation and the non-linear Reynolds Equation.McCarthy et al. [6] studied cantilever microswitches using a transient finite difference method approximating a parabolic pressure distribution along the length and non variance along the width and obtained good agreement with experimental measurements.Younis et al. [7] studied the effect of squeeze film damping for an electrically actuated micro-plate, using the compressible Reynolds Equation.They used perturbation methods to derive analytical expressions for pressure distributions in terms of the structural mode shape.Pandey et al. [8] studied the effect of flexural mode shapes on the squeeze film offered stiffness and damping for a cantilever resonator, they used Green's function to solve the linearized compressible Reynold's Equation and used the modal projection method available in ANSYS to solve the coupled fluid structure problem for several flexural modes of vibration.The analytical and numerical values of damping obtained were in good agreement with experimental results.Li et al. [8] accounted for the static bias deflection for a fixed-fixed micro-beam and a cantilever under electrostatic actuation.They assumed a parabolic function for the pressure along the beam width and a cosine series along the beam length, and solved the coupled Reynolds Equation and the Euler Bernoulli beam Equation.Hannot et al. [9] presented an approach to solve the coupled elasticity Equation and Reynolds's Equation for modeling a capacitive micro-switch.They employed a non-linear Newmark time integration scheme for the mechanical Equations and a trapezoidal rule for the fluid Equations.The above mentioned models attempt to solve the coupled problem, though not in a single step.The geometry modeled is also limited to 1D beam type structures.
These methods, though accurate, are cumbersome and involve iterative or staggered approaches.We propose a single step methodology to solve the coupled fluid-structure squeeze film problem.We use the 3D elasticcity Equation, thus not restricting ourselves to any particular geometry, and couple it with the 2D Reynolds Equation for squeeze film.A single step "monolith" approach [10] is presented to solve the coupled problem.The numerical results show good agreement with published experimental data and existing analytical solutions.

Numerical Modeling
The problem at hand involves solving for pressure on the vibrating plate due to the squeezed film, taking into account the elasticity of the plate.Thus the problem involves solution of the Reynolds Equation for the fluid domain coupled with the 3D elasticity Equation for the structural domain.In our finite element model, the air gap is treated as a 2D layer and the structural domain is modeled in three dimensions.The element used for modeling the structural domain is 3D, 27 node brick element.The "wet" face of the 3D element is treated as the fluid domain.Thus the relevant integrals for the fluid domain are evaluated over the corresponding 9 noded "wet" surface of the 3D, 27 node brick element.

Variational Formulation for the Fluid Domain
The linearized Reynolds Equation is given as follows [2], where μ eff is the effective viscosity, 0 h is the initial air gap, a P is the ambient air pressure, P is the fluid pressure (perturbed about a P ) and H is the air gap (perturbed about 0 h ).The last term on the right hand side of the above Equation couples the structure and the fluid domain.Substituting for harmonic solution in Equation ( 1) and considering p  as variation of p  in a weighted integral sense we have,

FEM Formulation for the Fluid Domain
For the FEM formulation we use 9 noded quadrilateral elements for 2D fluid domain.Pressure, its variation and z u  at any point are obtained from interpolation of the corresponding nodal values using the following relationships.
  Pressure gradient can be expressed as .

B B N N p N N u
(10)

FEM Formulation for the Structure
We have the variational statement for dynamic structural problem without anybody force as: where  is density, u is displacement,  u is its variation,   τ u is stress, t is prescribed traction over the boundary t  and    ε u is given by with the summation convention over repeated indices.
For coupled squeeze film damping problem with structural interaction, the wet surface (the surface which constitutes the 2D fluid domain) is subjected to prescribed traction ˆ, p    t n then Equation ( 11) can be written as where the last term on the left hand side signifies coupling effect of the fluid over the structural domain.The standard FEM discretization for displacement and other quantities are Substituting above relations in Equation ( 12) and using arbitrariness of  u we have the discretized Equation for the structural domain as

Coupled FEM Formulation
For the coupled problem at hand the 2D fluid domain corresponds to the "wet" surface of the structural domain.
Thus coupling the fluid and structure domains we have (combining Equations ( 10) and ( 13)) We have used 27 noded brick element for the structural domain, whose wet surface (consisting of 9 noded square face) is modelled as the 2D fluid domain.

Results and Discussion
For validation of our FEM results we have compared our numerical results with analytical solutions reported by Siddartha et al. [11,12], as well as experimental results from work by Pandey et al. [8].

Modeling a Varying Flow Boundary Elastic Microplate
Siddartha et al. [11,12] studied the effect of varying flow boundary conditions on the squeeze film stiffness and damping for an all sides clamped micro-plate.The plate is considered to vibrate in its fundamental mode which imparts the flexibility effect.Analytical expressions have been derived for stiffness and damping forces on the plate (clamped at all sides) due to the trapped air film (subjected to different flow boundaries).We use 4 × 4 mesh for FEM modeling of the plate.We design two representative flow boundary conditions with our numerical scheme, namely the all four sides open ("OOOO") and the two opposite sides closed ("OCOC") cases.The plate is subjected to harmonic displacement boundary condition corresponding to its first mode shape.The resulting pressure distribution is integrated over the wet surface to get the force on the moving plate.The squeeze film spring (F s ) and damping (F d ) forces are obtained from the real and imaginary component of the resultant force respectively.The forces are non dimensionalised (see [11]) and plotted against a non dimensional parameter , directly related to the frequency of harmonic excitation as follows, , where μ eff is the effective viscosity [13], ω is the harmonic excitation frequency, L is the plate side dimension, p 0 is the ambient pressure and h 0 is the initial air gap. Figure 2 shows the plots for F s and F d for the "OOOO" case and Figure 3 shows the same for the "OCOC" case.We see from these plots that the numerical stiffness and damping forces are in close agreement with the analytical results for both the flow boundary cases studied.We also note that the two methods show good agreement at both high and low  values (thus high and low frequencies).The deviation between the numerical and analytical results have been found to be less than 2% for both the flow boundary conditions studied here.

Modeling a Cantilever
In order to compare our results with experimental data we have modeled a cantilever beam as per dimensions  mentioned in the work of Pandey et al. [8].We compare the numerically obtained Quality factors (Q) for the first three modes of vibration as well as the effect of aspect ratio on the quality factor of the beam for the first mode of vibration.The beam modeled is of length 350 μm, width 22 μm, thickness 4 μm, with an initial air gap of 1.4 μm.The beam material is polysilicon with density 2330 Kg/m 3 , Young's Modulus 160 GPa and Poisson's ratio 0.22.Air is considered to be the fluid medium with the relevant property values (under standard temperature and pressure ) as follows: density ρ air = 1.2 Kg/m 3 , dynamic viscosity µ air = 1.8 × 10 −5 N.s/m 2 , and ambient pressure p a = 1.013 × 10 5 Pa.In our simulations we have subjected the cantilever to a sinusoidal voltage of 1.5 V, which is well below the pull-in voltage of 6 V for the given cantilever.The input voltage has been applied as an electrostatic pressure load of magnitude 5.08 N/m 2 to our FEM model for the cantilever beam, considering a parallel plate capacitor with small displacement approximation [14].The beam tip velocities have been obtained from the simulations and normalized with respect to the applied voltage and plotted against frequency.The frequency response so obtained is shown in Figure 4.The plot shows three distinct peaks corresponding to the first three modes, and is in close agreement with a similar plot reported by Pandey et al. [8].Q factors for different modes are obtained using half power method from the frequency response plot.Convergence study of the Q factor (Table 1) has been done using three levels of mesh refinement considering a very fine mesh (100 × 6 × 4) result as our benchmark.Q factors for sufficiently fine mesh (40 × 5 ×4 ) are compared with published results from Pandey et al. [8] in Table 2.We see that the data from the numerical simulations are in good agreement with published experimental and numerical results.We further studied beams with varying aspect ratios.We considered beams with lengths varying from 150 µm to 350 µm, having a constant width of 22 µm and thickness 4 µm.Only the first mode of vibration is considered in this study.Comparative values of Q factors for the different beams (40 × 5 × 4 mesh) are presented in Table 3.Our simulation results show a deviation of less than 10% from the reported experimental data.

Conclusion
We have discussed an FEM formulation to solve the coupled fluid-structure squeeze film problem without resorting to iterative solutions.Our results show good agreement with experimental data available from the literature.Our numerical scheme is seen to give good results for varying aspect ratio structures with dimensions below 100 μm.The proposed scheme can be further used as a design tool for modeling and simulation of the dynamic response of vibratory MEMS devices such as capacitive microphones, RF (Radio Frequency) MEMS switches, etc., for which accurate knowledge of the Q  factor is critical for design.With this methodology one can directly couple the elasticity effect of the structure with the fluid flow and need not limit oneself to 1D geometries.

For
on the closed borders.After doing integration by parts and implementing the above boundary conditions, we get the governing Equation for fluid domain as,

Figure 2 .
Figure 2. Spring and damping forces vs sigma for "OOOO" configuration.