Methods for Solving Gas Damping Problems in Perforated Microstructures Using a 2D Finite-Element Solver

We present a straightforward method to solve gas damping problems for perforated structures in two dimensions (2D) utilising a Perforation Profile Reynolds (PPR) solver. The PPR equation is an extended Reynolds equation that includes additional terms modelling the leakage flow through the perforations, and variable diffusivity and compressibility profiles. The solution method consists of two phases: 1) determination of the specific admittance profile and relative diffusivity (and relative compressibility) profiles due to the perforation, and 2) solution of the PPR equation with a FEM solver in 2D. Rarefied gas corrections in the slip-flow region are also included. Analytic profiles for circular and square holes with slip conditions are presented in the paper. To verify the method, square perforated dampers with 16–64 holes were simulated with a three-dimensional (3D) Navier-Stokes solver, a homogenised extended Reynolds solver, and a 2D PPR solver. Cases for both translational (in normal to the surfaces) and torsional motion were simulated. The presented method extends the region of accurate simulation of perforated structures to cases where the homogenisation method is inaccurate and the full 3D Navier-Stokes simulation is too time-consuming.


Introduction
Perforations are used in oscillating microelectromechanical (MEM) sensors and actuators to control damping due to gas for desired frequency-domain operation. Practical MEM devices using perforations include microphones [1], capacitive switches [2,3], and resonators [4].
The perforations reduce damping in a squeeze-film structure considerably: the characteristic dimensions are reduced from the outer dimensions of the surface to the separation of the perforations. Also, the cut-off frequency (when the spring and damping forces are equal) rises considerably due to the perforations. Often this cut-off frequency in perforated structures is much higher than the operating frequencies of interest, which simplifies the design of such devices since only the damping forces need to be considered in the device design.
3D flow simulations are needed in solving perforation problems generally. However, accurate 3D simulation of the Navier-Stokes (N-S) equation is, in practice, impossible due to the huge number of elements needed, especially in cases where the number of perforations is relatively large. The number of perforations may vary from a few holes to a grid of thousands of holes. Moreover, the N-S equations are limited to the region of a slightly rarefied gas (slip-flow region). Accurate modelling of these structures requires the consideration of rarefied gas flow in tiny flow channels and around the structure.
Different methods that reduce the 3D flow problem into 2D are available, depending on the boundary conditions, number of perforations, the relative hole size s 0 /h (s 0 is the hole diameter and h is the height of the air gap) and the perforation ratio (ratio of the perforated area and the surface area). If the s 0 /h-ratio is large, the flow resistance of the perforation is insignificant (unless the perforations are relatively long). In the case of a few relatively large perforations, the problem reduces to a 2D problem that can be solved from the Reynolds equation by considering the holes as additional boundary conditions [5,6]. Open edge corrections [7] might need to be considered in this case. On the other hand, if the number of holes is large, and the perforation ratio is not small, the horizontal flow from the damper borders is insignificant (also caused by closed damper borders), and the pressure patterns around each perforation are identical. In this case, the damping can be determined solely from a single perforation and its surroundings, called a perforation cell. A simple analytic model for such a cell by Skvor [8] has been used in modelling perforations in a microphone's backplate. When the s 0 /h-ratio cannot be considered large, the flow resistance of the perforations must be included in the model. Homentcovschi and Miles [9], have derived such an analytic model considering the compressibility and inertial effects. Perforation cell models over a wide range of perforation ratios and hole lengths well below the cut-off frequency are given in [10,11].
If the s 0 /h ratio and the perforation ratio are not large, such that the flow from the damper borders is considerable, a more complicated model is required. In the case of a regular grid of perforations, the flow resistance of the perforation cell can be solved and homogenised over the damper surface. This flow problem can be modelled with the extended Reynolds equation that contains an additional pressure leakage term [12][13][14][15][16]. This homogenisation method is usable, but it has its limitations, especially if the number of perforations is small. A mixed-level approach has also been published [17], and a hybrid method combining N-S and Reynolds equations has been presented by da Silva et al. [18]. Flow resistances of perforations with different shapes were studied by Homentcovschi and Miles [19]. A compact analytic model for rectangular dampers having an uniform perforation was given in [10]. The model was verified with FEM simulations at frequencies well below the cut-off frequency.
The models referenced above assume motion normal to the surfaces, and the squeezing force in the air gap makes the gas flow. A different problem where the perforated surface moves tangentially to the surfaces has been discussed in [20]. We have presented a Perforation Profile Reynolds (PPR) method [21] for damping problems where the homogenisation method is not accurate. These are characterized by, e.g., a moderate number of perforations, complex shaped perforations, and non-uniform distribution and sizes of perforations.
In this paper, we present improvements to this method. Flow admittance profiles that vary across the perforations are used instead of constant profiles. Also, an improved perforation cell model [10] has been used instead of a simple theoretical long-channel model with theoretical end corrections. Besides translational motion, the verifications include also torsional motion. This paper also shows how the perforation cell model [10] can be utilised in a homogenisation method with a FEM solver.
The method was verified with 3D FEM simulations using a N-S solver. Simulations of a parallelplate damper with 16 to 64 square perforations were performed for various perforation ratios. Only incompressible flow is considered here, since time-harmonic reference N-S simulation with sufficient accuracy is not realistic. Figure 1 shows the general geometry of a perforated gas damper. A perforated surface above the ground surface moves normal to the surfaces (in the direction of the z-axis). Gas lubrication and damping problems in small air gaps are traditionally modelled with the Reynolds equation [22]. The modified and linearised Reynolds equation in the time-harmonic form is

PPR equation
where h is the static air gap height, η is the viscosity coefficient, Q ch is the relative flow rate, P A is the ambient pressure, v z (x, y)e jωt is the surface velocity in the z-direction and p(x, y)e jωt is the pressure variation to be solved from the equation. This form allows the static displacement h(x, y) to be a function of x and y. In this case, the spatial dependency of Q ch (x, y) must be considered, too. Isothermal conditions have been assumed.
We extend the Reynolds equation by adding a term that models an additional flow component through the perforations similarly as in [12]. Here, two additional coefficients are included to model the effect of the perforations on the diffusivity and compressibility.
where D h (x, y), C h (x, y), and Y h (x, y) are the relative diffusivity, relative compressibility, and perforation admittance profiles, respectively. Equation (2) is called here the Perforation Profile Reynolds (PPR) equation.
The perforation admittance Y h is specified as an inverse of the specific acoustic impedance Z S : where p h and v h are the gas pressure and velocity in the z-direction at the opening of the perforation. Equation (3) is valid also for the quantities averaged across the cross-section of the perforation, denoted as Y h , Z S , v h , and p h . The approximation for the relative flow rate due to the rarefied gas effects is [23] Q ch = 1 + 9.638K 1.159 n,ch , where K n,ch = λ/h is the Knudsen number of the air gap and λ is the mean free path of the gas (depends on P A ). The relative accuracy is valid for K n,ch < 880 and its maximum relative error is less than 5%. When the slip correction is considered, as in this paper, The slip correction is valid only for K n,ch < 0.1, whereas the approximation in Eq.( 4) is valid for much larger Knudsen numbers. The slip correction is used in this paper, since for us the validation with a 3D FEM solver is only possible in the slip-flow region.

PPR method
The PPR method consists of two phases. First, diffusivity, compressibility, and perforation admittance profiles in Eq. (2) are determined. Then, the PPR solver is used to solve Eq. (2). In PPR, the specific admittance profile Y h (x, y) varies spatially at the locations of the perforations, whereas in the homogenisation method Y h includes the flow resistance of the perforation cell, that is, the flow in the air gap and the perforations.
In [21], we used a constant average admittance profile for each perforation. This admittance was calculated from the flow resistance of the perforation. In this work, we consider a varying admittance profile across the perforation utilising the velocity profile in the perforation: where v h is the velocity profile and p h is the pressure. The diffusivity of the air gap also effectively changes at the perforation. For the flow passing in the air gap below the perforation, the frictional surface area will be reduced, effectively increasing the diffusivity coefficient. On the other hand, the flow profile is changed due to the perforations, effectively decreasing the diffusivity coefficient D h . The compressibility coefficient C h is zero at the perforations and unity elsewhere. Figure 2 illustrates the perforation profiles in the case of a one-dimensional damper.

Perforation Admittance Profiles
There are two alternate ways to determine the flow admittance profiles for perforations: 3D FEM (or other discretising methods) simulations and analytic (approximate) equations. Usually, the perforations are relatively short and the analytic model for long channels are not applicable. Here, flow resistances extracted from FEM simulations of the perforation cells are used [10,11].
The flow fields above the perforated surface are coupled, especially when the perforation ratio is large. In [10] an elongation model considering the coupling in circular perforations is derived from FEM simulations. The model assumes that the conditions for neighbouring perforations are identical. In principle, these conditions are met if the perforation ratio is large and the velocity distribution v z (x, y) of the damper surface is constant. However, in practice the coupling is negligible for small perforation ratios, and there are no abrupt changes in the damper surface velocity for neighbouring perforations.
Next, analytic expressions for two typical perforations, circular and square, are given (See Fig. 3).

Circular Cross-section
For the flow resistance of the perforation we utilise the perforation cell model [10]. Here, the flow resistances modelling the perforation only are needed, since the PPR solver calculates the flow in the air gap. The (average) specific admittance of a channel (of length h c ) with a circular cross-section (with a radius of r 0 ) in the slip-flow region is [10] where the flow resistances R IC , R C and R E are given in Appendix. R C is the mechanical resistance of a long channel with a circular cross-section. In the slip-flow region, the relative flow rate of a channel with a circular cross-section Q tb is where K n,tb is the Knudsen number of the capillary The normalised velocity profile (average velocity is 1) in a long circular capillary is a function of r only [24]: The specific admittance profile resulting from Eqs. (7) and (10) is

Square Cross-section
Mechanical flow resistance of a long square channel is [25,26] where Q sq is the flow rate coefficient and h c is the length of the channel. Next, the short-channel effects are added to R sq using the elongation model for a circular channel [10]. In that model, an effective radius r 0 is needed. It is specified such that the acoustic resistances (hydraulic resistances) of the square and circular long channels are equal. That is This results in an expression for r 0 : Since Q tb is a function of r 0 , the accurate solution of Eq. (14) needs to be iterated. The dependency of r 0 on the Knudsen number is small, and the simple approximation in Eq. (14) is accurate within 0.3% for K n,sl < 0.14.

Homogenisation method
In the homogenisation method, the flow admittane Y h caused by the perforations is assumed to be smoothly distributed over the whole surface. In the actual 2D FEM simulation, the holes are excluded from the simulated structure and their flow admittance is included in the extended Reynolds equation. In the case of uniform perforation, this flow admittance is constant over the whole surface. Damping models using the homogenization approach have been published, but not properly verified [12][13][14][15][16].
The PPR equation (2) is applicable also in the homogenisation method, but there is a considerable difference in the interpretation of term Y h . In addition to the flow admittance of the perforations itself, Y h should now include the flow component in the air gap.
In deriving the model for Y h , the concept of a "perforation cell" is applied here. Relative simple analytic models for a cylindrical and square perforation is derived in [10] and [11], respectively, based on compact models derived from FEM simulations. The mechanical resistance of each perforation cell consists of mechanical resistances of several flow regions [10]: where r x and r 0 are the outer and inner radii of the cylindrical perforation cell, respectively. R S is the flow resistance in the air gap and R C is the flow resistance of the perforation. All other resistances model the flow in the intermediate region under and above the perforation. Values for these resistances are given in Appendix.
Radius r x is selected such that the bottom areas of circular and square cells match. This results in In the case of uniform perforation and translational motion, the specific impedance spread all over the surface is where a eff = a + 1.3h(1 + 3.3K n,ch ) is the effective dimension of the surface.

Verification
Several perforated dampers in perpendicular and rotating motion are simulated both with the 3D N-S solver with slip conditions and with the time-harmonic PPR solver. Both solvers are implemented in the multiphysical simulation software Elmer [27].
The surface is a square, and it consists of N square holes. The distance between the edges of two adjacent holes is s 1 , as is the distance from the edges of the holes to the edges of the damper s 2 = s 1 . The perforation pitch is s x = s 0 + s 1 . A comparison is made well below the cut-off frequency, that is, the compressibility of the gas is ignored (C h = 0). The simplifies the N-S simulations, since response to a constant velocity is calculated. The dimensions are summarised in Table 1 and one of the four structures is shown in Fig. 4.

3D simulations
First, full 3D N-S simulations of the structure are performed to generate reference data for the methods presented in this paper. The simulations can be performed relatively reliably, thanks to the small number of holes and the symmetry of the structure. Slip boundary conditions are used for the surfaces. The simulated gas volume is extended around the damper: free space around and above the damper are 6h and 2s 0 , respectively. A mesh of 450 000 elements is used. Both translational and torsional motion cases were simulated.
The number of elements (450 000) in the 3D N-S simulations was in practice the maximum that could be performed considering the computing time and the memory consumption. The computer used was Sun Fire 25K. Simulating each (of 432) topology took between one and two hours and 6 Gb of memory.
To estimate the accuracy of the results, simulations with 250 000 elements were made, too. Comparative results are discussed later in this paper.

PPR solver
The PPR equation (2) is solved with a time-harmonic PPR solver. For each square hole, the specific admittance for a circular hole specified in Eq. (11) with an effective radius of r 0 is used. The PPR solver also includes the compressibility of the gas in the air gap. The compressibility of the perforations can also be included specifying a complex-valued specific admittance. This is essential in modelling, e.g., closed perforations. Compressibility profile is not needed here, since steady flow calculations are made.
We use D h = 1 since there is no model available for the diffusivity profile. The friction surface for the flow under the hole is reduced by half, and this could justify setting D h = 2. However, the changes in the flow profile will introduce additional losses that justify decreasing the diffusivity.
Considering the edge effects at the outer borders of the structure is essential for accurate results. Each border of the structure is extended by an amount 0.65h(1 + 3.3K n,ch ), as suggested in [7]. The borders of the perforations are not modified in this method, since the diffusivity profile at the perforation effectively builds an "edge effect": the flow admittance at the border of the perforation is small, forcing gas to flow horizontally in a viscous channel close to the border.
Both translational and torsional motion cases are simulated. Each of the simulations with the Reynolds solver with 20 000 elements took about 10 seconds.

PPR method
First, the structure is set into translational motion with a velocity of v z = 1 m/s. The mechanical resistance, or the damping coefficient is Z m = F/v z , where F is the force acting on the ground surface. A mesh of 20 000 elements is used.
The damping coefficients are compared with the ones from the N-S simulations as a function of the perforation ratio The maximum relative errors in the damping coefficient are shown in Table 2, and Fig. 5 shows the damping coefficients as a function of the perforation ratio q for air gaps of 1 µm and 2 µm. In Fig. 6, the simulated pressure profiles are compared. Next, the structure is set into torsional motion about the y-axis with a velocity of v z = 1 m/s at x = a/2. The torsional resistance, or the torsional damping coefficient is Z t = τ/ω z , where τ is the twisting moment acting on the ground surface and ω z = v z /(a/2) is angular velocity. A mesh of 20 000 elements is used. The maximum relative errors in the torsional motion are shown in Table 3. Figure 7 shows the torsional damping coefficients as a function of the perforation ratio for air gaps of 1 µm and 2 µm.

Homogenisation method
The extended Reynolds equation is solved for perpendicularly moving rectangular surfaces. The perforation cell admittance calculated from Eq. (17) is used in the simulation. The diffusivity coefficient D h = 1 is used. Each damper border is extended by 0.65h(1 + 3.3K n,ch ). The quarter of the rectangular surface is meshed having 5 000 elements. Table 4 shows the maximum relative errors of the homogenisation method.   Table 2 shows the errors of the PPR method compared with the 3D FEM simulations for translational motion. The maximum relative error in the damping force in translational motion is below 10%, except for the case of strong perforation (q ≥ 49%). The PPR model generally underestimates the damping. This is probably due to the fact that additional drag forces acting on the sidewalls and on the top surface are ignored in the model. These forces contribute most when the surface is the thickest (h c = 5 µm) and the surface area is the smallest (N = 16).

Comparison between methods
In the case of torsional motion, Table 3, the accuracy of the model is approximately the same as for translational motion for N = 64. The contribution of the drag forces can be seen in the results of smaller surfaces in Fig. 7. In the torsional motion the drag forces acting on the sidewalls are relatively stronger. For small q, the change in the height of the surface changes the damping force for N = 16 in the N-S simulations, see Fig. 7b). The surface height clearly changes the damping coefficient (N-S simulations) even for very small perforation ratios. This change cannot be explained by the change in the perforation length, since the damping due to perforations is negligible for very small perforation ratios.
To analyse the importance of the varying pressure profile at the perforations, a set of PPR simulations was performed using a constant specific admittance Y ci of the perforations. The results are summarised in Table 5. Comparison with results in Table 2 reveals that errors of both models are approximately the same for small perforation ratios, but the error of the constant specific admittance model is considerably larger if q > 16 %.
The maximum relative errors of the homogenisation method (Table 4) are also approximately below 10 % for N ≥ 16 and q < 49 %, but the PPR results are much better for small perforation ratios. However, the homogenisation method gives better results when N ≥ 16 for q ≥ 50%. This was an expected result. For an inhomogeneous perforation, where the distribution and sizes of the perforations vary, the determination of the perforation cell might be difficult. In such a case, the generation of the variable perforation admittance profile is not a trivial task. On the other hand, if the number of perforations is large (hundreds or thousands of holes) the element sizes in the PPR method will grow intolerably. In this case, the homogenisation method is superior since the required element count is independent of the number of perforations.

Rarefied gas effects
The verification in this paper is limited to the slip-flow region just to make the comparison with 3D solutions possible. The PPR method itself is valid for arbitrary gas rarefaction as long as the model for Y h is valid. The comparison is made in the region where the validity of the slip conditions is questionable (K n,sl < 0.14, K n,ch < 0.14). In spite of the small additional error in the results, the comparison is justified since slip flow models are used in both cases. The consideration of the slip correction is essential for accurate results for the micromechanical structures simulated. If the rarefied gas effects are ignored (K n,sl = K n,ch = 0), the error in the damping force would increase by about 40 % and 15 % for air gaps of 1 µm and 2 µm, respectively.

Maximum frequency of the model
In squeeze-film dampers, the compressibility of the gas will change the damping coefficient at higher frequencies. The cut-off frequency specifies a frequency when the viscous and compressibility (spring-) Table 5. Maximum relative errors in the translational motion damping coefficient Z m simulated with the PPR method. Average specific admittances are used for the perforation and D h = 1. All other simulation parameters are the same as in the simulations in Table 2. forces are equal. For perforated dampers, the small-frequency assumption made in this paper is generally a valid, since the perforation increases the cut-off frequency considerably. The cut-off frequency due to compressibility can be estimated by calculating the ratio between the viscous and compressibility forces in a perforation cell [10]. This is quite straightforward after assuming that the perforation ratio is high enough that all gas flows through perforations (not from the damper borders). However, the cut-off frequency might be so high that the inertial effects should also be considered. In estimating the validity the complex flow resistance due to inertia should be considered. In [10] the maximum frequency for a perforation cell is estimated.
There is no limitation in the method that prevents compressibility from being accounted for. In this case, a complex-valued specific admittance profile and compressibility profile C h are both needed for the PPR solver. For verifying this case a 3D linearised time-harmonic N-S solver may be used.

Accuracy of 3D N-S simulations
The 3D simulations were very challenging. A lot of computing resources were used to have the data, still the results are not perfect. To estimate the error in the simulations with 450 000 elements, a comparison is made to simulations with 250 000 elements. The results are shown in Table 6.
When the perforation ratio is large, the relative difference is the largest. This indicates that the selection of the meshing was not very good for large perforations.
The comparison shows that damping coefficients with 250 000 elements are generally larger than with 450 000 elements. It is then expected that the damping coefficients become smaller when the element count is increased. When considering such fictitious simulations and studying the maximum relative errors of the PPR method in Table 2, the magnitudes of the relative errors become generally smaller.

Conclusions
A straightforward method solving perforation problems with a 2D PPR solver was presented. Damping coefficients of a large number of perforated dampers were solved with the PPR method and they were compared with 3D N-S simulation results with very good agreement. Compared to direct N-S simulations, the PPR method offers orders of magnitudes faster simulation times and less memory consumption. The maximum relative errors for translational motion for N ≥ 16 were about 10 % and 20 % for perforation ratios 1%. . . 36% and 1%. . . 81%, respectively.
Perforations were modelled with their specific admittance profile. Both, constant and spatially varying admittance profiles at the perforations were used. It was shown that the errors were smaller when a varying admittance profile was used. Uniform diffusivity profile (D h = 1) was assumed. A further study is needed to resolve if a non-trivial profile could be extracted form FEM simulations, and if the utilization of this profile would improve the accuracy of the model.
The PPR method was also compared with the homogenisation method. The results were almost as good (maximum relative error about 10%) for N ≥ 16 for small perforation ratios (q < 36 %).
The square perforations were approximated with a model for a cylindrical perforation cell. It is expected that a model for a square cell [11] could improve the accuracy of PPR and homogenisation methods, especially at high perforation ratios.
The fringe flow effects due to the sidewalls and the upper surface clearly contributed the damping coefficient for small s 0 /h ratios. This was seen in the results for dampers with a small number of holes, since in these cases the surface area was also relatively small.
Incompressible flow was assumed. For a squeeze-film damper this sounds like a very limiting assumption. But due to the perforations, the cut-off frequency is significantly higher than in the non-perforated case. Thus the valid frequency range of the incompressible damper model is relatively wide. Further increase of the frequency range is not trivial, since, besides the compressibility effects, the inertia of the gas should be considered. Also, the thermal aspects should be considered since isothermal assumptions will not be sufficient due to the high cut-off frequencies.
The method needs to be verified in cases where the gas compressibility and inertia are considered. However, this is not easy, since the time-harmonic analysis of the 3D structure with a sufficient accuracy will be very challenging. Extraction and usage of non-uniform perforation profiles and coupling between perforations will be studied in the future.  Figure 8 shows the topology of the cylindrical perforation cell, and Fig. 9 shows the lumped mechanical resistances used in modelling the flow resistance of the cell. The flow resistance R p of the perforation cell consists of lumped flow resistances and their effective elongations. The equations have been derived partly analytically and partly by fitting the model to FEM simulations by varying the coefficients in heuristic equations [10].

Appendix
The mechanical resistance of the perforation cell is [11] R p = R S + R IS + R IB + r 4 where r x and r 0 are the outer and inner radii of the cylindrical perforation cell, respectively. The lumped resistances are