Dynamics of a Micro Electrical Mechanical System Subject to Thermoelastic and Squeeze-Film Damping

The static and the dynamic behaviours of a MEMS subjected to thermoelastic and squeeze film effects have been investigated. Major attention is initially devoted to the modeling of this multi-physics system, including mechanical, electrical, thermoelastic and fluid dynamics behaviours, and their couplings. Then, static solutions have been obtained numerically by a finite-difference technique. Finally, an analysis of the vibrations is performed, and the complex natural frequencies are determined. The real part is related to the oscillatory behaviour, while the imaginary part provides the overall damping of the system.


Introduction
The analysis of the linear and nonlinear dynamics of Micro-and Nano-Electrical-Mechanical-Systems (commonly known as NEMS and MEMS, respectively) have been largely investigated in the recent past [1] due to the increasing interest toward this particular class of dynamical systems.
The practical importance arises from a variety of applications in challenging areas of research and development, where using MEMS enables very smart solutions unavailable with other technologies.Some examples of MEMS are accelerometers [2], inertial sensors [3], chemical sensors [4], filters and oscillator [5].
From a theoretical point of view, on the other hand, the interest in focused on modeling [6][7][8][9] and on the analysis of the response [10,11].In the latter case the goal is often that of accurately determining [12] and controlling [13] the static and the dynamic pull-in threshold, which could or could not be an unwanted phenomenon depending on the application at hand, and in fully understanding the overall dynamical behaviour in order to improve the performances of the system [1].
One of the crucial points of the modeling activity is related with the multi-physics nature of the system, which includes mechanical behaviour and electrical actuation, which are the basic "ingredients" of any MEMS, and possibly thermoelastic [14][15][16] and squeeze film (fluid dynamics) [17,18] effects, which play a key role in certain applications, such as, for example, in MEMS operating in non vacuum environments [19].
Thermoelastic and squeeze film effects are the main sources of energy dissipation, and thus they significantly affect the quality factor (which is inversely proportional to the overall damping) and the performances of MEMS.Thus, when present in real cases, they need an accurate modeling in order to properly predict the effective behaviour.
Although they have been repeatedly analyzed [1], to the best of our knowledge they have been never considered jointly, both in modeling and in simulations, a fact that constitutes the main goal of this paper, which is a first report of a more general work which is still in progress.

The model
We consider the rectangular microbeam illustrated in figure 1, which is supposed to be fixed at both ends.In between the beam and the substrate there is a fluid, whose dynamics is influences by the motion of the beam.This section is devoted to the determining to the coupled partial differential equations governing the planar dynamical behaviour.Since there are 4 fields (mechanical, electrical, thermoelastic and fluid dynamics) involved in the problem, we divide the section in 4 subsections.

Mechanical behaviour of the beam
From a purely mechanical point of view, the beam is considered within the Euler-Bernoulli approximation, i.e. unshearable and with cross-section remaining orthogonal to the axis of the beam.
The governing equation is then [20] ̂ where is the transversal displacement, is the mass density, b and d are the cross section dimensions (see figure 1), EI is the bending stiffness, ̂ is the mechanical damping coefficient (measuring the internal damping of the beam), is the axial stiffness, is the length of the beam, is the axial load.Furthermore, we have that is the distributed moment due to temperature, and are the transverse loads per unit length due to the electrical and fluid dynamics effects.These three latter terms represent the coupling of the mechanical behavior with the others three physical fields.
In (1) the axial inertia has been neglected, so that the axial displacement has been condensed in the classical way, the material is linear elastic, and the geometrical nonlinearity is considered only in computing the elongation.It is expanded up to the third order, and it is represented by the integral term between brackets.
The boundary conditions are , at and . (2)

Electrical behaviour
The system constituted by the beam, which is flat in the unstressed configuration, and the lower electrode (substrate), which is placed directly underneath the microbeam at a rest distance , is a capacitor.The time varying gap is assumed to be smaller than the beam length L. This entails the electrical field to be approximately constant along the beam axis.With this hypothesis, the field equation can be solved in closed form, leading to where , is the dielectric constant in the free space, is the relative permittivity of the gap space medium with respect to the free space, is the voltage difference between the beam and the electrode, and being the electrostatic and the electrodynamic voltages.
The voltage is assumed to be known, so that the electric field does not introduce a new equation in the problem.Eq. ( 3) is strongly nonlinear with respect to the unknown , this being at the base of the most important mathematical difficulties arising with the MEMS dynamics.
The pull-in condition corresponds to an infinite force acting on the microbeam, which is clearly a mathematical limit of a physical situation in which the beam gets in contact with the substrate.In some cases this is an unwanted situation (short circuits), which could even destroy the MEMS, while in other cases is it useful (e.g. in switches).
Since the length of the beam is assumed to be large enough, we neglect the fringing effects of the electric field at the edges of a plate capacitor.

Thermoelastic behaviour of the beam
In some types of MEMS the loss of energy due to thermoelastic effect could be important.The phenomenon is strictly related to the frequency of the vibrations.When the period of the oscillations is much larger or much smaller than the thermal relaxation time, we are in isothermal or in adiabatic frameworks, respectively, and the thermal losses are often negligible.
For comparable periods, thermoelastic damping is a fundamental dissipation mechanism that is inherent to the system and cannot be eliminated, and thus must be carefully modeled in order to detect its effects.
The constitutive law for an isotropic thermoelastic beam is , where is the Young modulus, is the axial elongation of the axis of the beam, the curvature, is the thermal stress coefficient, is the linear thermalexpansion coefficient, is the difference between the actual, , and the reference, , temperatures.
From (4) it is easy to compute We need the evolution equation for , since it is neither constant nor known.The thermodynamics laws and ( is the work, is the internal energy, is the heat, and is the entropy) provide ( and are the stress and strain tensors, respectively) We have ( ), ( ) and ̇ , where internal heat sources are disregarded, is the local heat flux density, and the dot represents derivative with respect to time.
According to the Lord-Shulman [21] or Maxwell-Cattaneo [22] theories, in the isotropic case and are related by where is the relaxation time and is the thermal conductivity.Clearly ( 7) is an extension of the classical Fourier law, which is obtained for .After some manipulations we get from the previous equations where is the specific heat capacity.Assuming small we get In the beam we have only , so that ( 9) becomes To further simplify (10), we assume that, due to the small thickness of the beam, the temperature profile can be approximated as which fulfills the boundary conditions at .Inserting (11) in (10), multiplying by , integrating in from to and remembering the definition (5) of we get where ∫ .The previous expression is the governing equation in the thermal bending unknown .

Dynamics of the fluid between the beam and the substrate
In between the two plates of the capacitor there is a fluid, which actually is a gas, commonly the air.When the beam moves down toward the substrate tends to "squeeze" out the fluid, while when it moves up it tends to "suck" in the fluid from the lateral sides of the capacitor.
To model this phenomenon, we initially assume that along the gap ( -direction) the velocity is parabolic, and that the fluid always remains in contact with the beam during the motion.
The pressure of the fluid is governed by the Reynolds equation [17] where ̅ is the effective viscosity and , ̅ , ( are the time varying gap and the time fluctuating pressure, being the initial static pressure, which could be the atmospheric one in some cases.In (1) it is ̅ .Inserting ( 14) in (13) and taking only the linear terms in ̅ and we obtain If we further neglect the dependence of ̅ with respect to , since the length in the direction is much smaller than the length in the direction, we get which is the governing equation for ̅ .

Governing equations
Summarizing the results of the previous subsections, the system of three partial differential equations in the three unknowns , (which becomes in dimensionless form) and ̅ is, in dimensionless form, The associated boundary conditions at and are instead , , and

Results
In this section the static deflection and the free linear oscillations are considered.The dynamic voltage, which constitutes the excitation, is, therefore, vanishing, .

Solution strategy
To start with the steady-state case, we neglect the time derivatives in (17), obtaining Solving the second and the third equations together with the boundary conditions ( 18) provides , , and then is the given the static pressure of the fluid between the beam and the substrate, and practically corresponds to the difference of pressure below and above the beam.In the following we assume that this difference is zero, .Equation ( 20) has been solved numerically by a finite-differences scheme and by a Newton-Raphson selfmade algorithm.
As a second step we compute the linear oscillations around the nonlinear static configuration .With this in mind, we look for a time-varying solution in the form where is the complex frequency of oscillation.The real part is related to the period of oscillation, while the imaginary part is related to the overall damping.Inserting (21) in (17), neglecting the nonlinear terms, and equating to zero the coefficients of and of we obtain the systems of 6 homogenous ordinary differential equations: where prime means derivative with respect to .The solution of ( 22) is too much involved to be determined in closed form, and so it is sought in the form where is the n-th eigenfunction of the problem with boundary conditions .
Equation ( 24) is clearly a simplification of ( 22) 1,2 , and it is chosen because the solution can be written in closed form and, being self-adjoint, the are orthogonal.By applying the Galerkin method we arrive at a homogenous linear algebraic system in the unknowns , , , , , .Equating to zero the determinant of the associated matrix provides the characteristic equation whose solutions are the frequencies .The localization of the solution is quite difficult, and it has been obtained by using a corollary of the residues theorem.

A numerical example
In this section we consider a numerical example referring to the physical parameters of Table 1.
The dimensionless parameters appearing in (17) follows from those of Table 1 and are reported in Table 2.
In the following simulations we keep fixed all the parameters of Table 2, apart from the damping and the coupling parameters and which varies from 0 to the values reported in  Since the precision of the solution obviously depends on the number of the considered eigenfunctions in (23), a convergence analysis is illustrated in figure 2, where it is shown that with we obtain an adequate accuracy in the considered range.
The first five complex frequencies for increasing are reported in figure 3. The left ending points of each path are given in (27), while the right ending points in (28).The five paths are regular and qualitatively similar, although that of is more curved than the others.When the coupling increases, the imaginary part increases as well, so that the overall damping increases and the quality factor decrease, as expected.With the present analysis we are able to quantitatively determine the reduction of the performances of the MEMS.
The real part increases as well, meaning that the thermoelastic and fluid dynamics effects strongly influence also the oscillatory behavior.This is especially true for low frequencies, while for high frequencies the paths tend to become vertical, i.e. the lower periods weakly depend on the damping.
It is worth to remark that, for fixed values of the parameters, the overall damping increases with the order of the frequency.
From figure 3 we see that the paths of and get closer and closer, suggesting that for larger parameters an internal resonance might be possible.

Conclusions and further developments
The governing equations of a generic MEMS device subjected to thermoelastic and squeeze film damping have been obtained.The multi-physics nature of the CSNDD 2012 04004-p.5 problem has been stressed, and the coupling terms analyzed.The (complex) frequencies of the damped free oscillations around a static nonlinear configuration have been determined by numerically solving the static nonlinear problem, and by the Galerkin method in the dynamic part.
The results presented in this paper are preliminary, since they are part of a work which is still in progress.
In addition to performing further simulations, for example by doing a systematic parametric analysis, by considering the influence of the nonlinear static deformation, etc., various other developments are worthy.Among them we think important a proper account of: 1) the thermoelastic behavior not only of the beam, but also of the fluid between the plates of the capacitor; 2) full three dimensionality of the problem, using for example the Kirchhoff-Love plate theory; 3) the nonlinear effects around resonance, by determining the amplitude dependent nonlinear frequencies and, in the forced case, the bending of the frequency response curve.
These are left for future works.

Fig. 2 .
Fig. 2. The convergence analysis for the first (top panel) and the second (bottom panel) frequency for increasing and for increasing number of considered eigenfunctions.

Fig. 3 .
Fig. 3.The first five complex frequencies for increasing values of .

Table 1 .
The physical parameters.