PDF
Squeeze-Film Gas Damping of a Vibrating Disc
Introduction
This benchmark model computes the damping force acting on a vibrating disc. The disc is in close proximity to a stationary surface and the damping results from the squeezing of a thin film of gas between the two surfaces. The squeezing action forces out the gas from between the two plates resulting in a damping force that acts to prevent mechanical contact between the two surfaces. The opposite effect takes place when the surfaces move away from each other as gas is drawn back into the bearing.
This model examines the effect of the periodic motion of the disc on the flow developed, including the pressure in the gas and the resulting damping forces. Small amplitude motion is analyzed using a linear frequency domain simulation. A nonlinear transient analysis is performed for small to large amplitude motion. The calculated film pressure and load carrying capacity are compared with analytical results.
Model Definition
The Thin-Film Flow, Edge interface is used to model the gas film on a flat circular plate. The model is 1D axisymmetric since the film pressure only varies radially. When Thin-Film Flow is assigned to a boundary, this boundary represents a reference surface in the physical device. In practice a small gap exists at the boundary and two impermeable structures, the wall and the base, are located either side of it. Figure 1 shows the configuration of the base and the wall in an arbitrary problem, and defines a number of terms used in the interface.
In this example, the model geometry is 1D axisymmetric and consists of a single line, with length set to the radius of the circular disc. The line is located at the origin and aligned with the r-axis. The base is coincident with the reference surface. A pressure is generated in the bearing by a periodic velocity of the wall in a direction normal to the wall.
Figure 1: An example illustrating a typical configuration for thin-film flow.
For non-slip boundary conditions at the wall and the base, the modified Reynolds equation takes the following form for a general frequency domain problem:
(1)
where ρ is the fluid density, μ is its viscosity, h0 is the mean film height, pf is the pressure developed as a result of the flow (this is the dependent variable in COMSOL) and ptot is the total pressure (ptot = pA + pf, where pA is the ambient pressure). Other terms are defined in Figure 1. A reference surface with normal nref sits in a narrow gap between a wall and base. In COMSOL the vector nref points into the base and out of the wall. The wall moves with a displacement from its initial position with displacement uwall and velocity vwall. Similarly the base moves from its initial position with displacement ubase and velocity vbase. The compression of the film results in an excess pressure, pf, above the ambient pressure, pA, and a gas velocity in the gap. At a point on the reference surface the average value of the film velocity along a line perpendicular to the surface is given by the in plane vector vave. The motion of the gas results in forces on the wall (Fwall) and the base (Fbase). The height of the wall above the reference surface is hw whilst the base is a distance hb below the reference surface. The total size of the gap is h=hw+hb. At a given point in time hw=hw1–nref · nwall and hb=hb1–nref · nwall where hw1 and hb1 are the initial heights of the wall and base, respectively.
Note that the frequency formulation assumes a small amplitude first order harmonic variation of film pressure, film height and wall velocity at the frequency of interest. The boundary conditions for this model are vanishing pressure due to the flow (pf = 0) at r = r0, where r0 is the disc radius and symmetry/zero pressure gradient (dpf/dr = 0) at r = 0.
For the case of a 1D axisymmetric problem Equation 1 can be greatly simplified to derive a simple closed form analytical solution. Note that these simplifications are not made in the simulation itself and consequently there are slight deviations from the analytic results in the model — these cannot be seen in the plots shown here and are not significant (strictly speaking the simulation is more accurate than the analytic results since no assumptions are made). The motion of the disc is in the vertical direction only, and the gap size is uniform across the disc, so a number of the terms in Equation 1 are zero. In this example, the term iωhpf is quite small compared to other terms and can be neglected for the purpose of deriving an analytical solution. The ambient pressure is 1 atmosphere and correspondingly pf«pA. so ptot≈pA. Making these assumptions the modified Reynold’s equation reduces to:
(2)
Where vw is the velocity of the wall in the z-direction. With the boundary conditions pf = 0 at r = r0, and dpf/dr = 0 at r = 0, Equation 2 can be solved for pf and is given by (see Ref. 1 for complete derivation):
(3)
The total analytical vertical load on the disc is then given by (see Ref. 1 for complete derivation):
A frequency domain analysis is appropriate for small amplitude periodic motion of the bearing wall. For large amplitude periodic motion (where the amplitude of the motion becomes comparable to the gap size), a transient analysis is necessary, to capture the nonlinearities in the model. Again approximations are required to derive an analytic result for the transient simulation (but are not necessary for the model).
For a transient model the modified Reynolds equation takes the form:
(4)
For periodic motion of the wall, the total film height h is h(t) = h0 + Δhsin(2πft), where Δh is the amplitude and f is the frequency of wall periodic motion. The wall velocity vw is then given by vw(t) = (2πfΔh)cos(2πft). Once again making the assumption that ptot ≈ pA and noting that a number of these terms are zero for vertical motion of a parallel disc:
Following the derivation of total analytical vertical load for the frequency domain analysis, the total vertical load on the disc for the transient model is given by
The model compares the analytical values of film pressure and total vertical load against the values computed by using the Thin-Film Flow, Edge interface. The results are found to be in agreement with the analytical solutions.
Results and Discussion
The values of radial film pressure in the gas film for the frequency domain analysis are plotted in Figure 2. As expected, film pressure magnitude is maximum at the center of the circular disc and drops off to zero where the gas film exits the bearing geometry. Figure 2 also compares the numerical radial film pressure values with analytical values calculated using Equation 3. The calculated results agree well with the analytical results. Figure 3 shows an arrow plot of the fluid load per unit area on the circular disc. These values correspond to the film pressure since the vertical load is significantly larger than the radial load. Figure 4 shows the variation of film height (gap) with respect to time for different values of the amplitude of harmonic film height. The corresponding variation of film pressure and total vertical load on the circular disc is shown in Figure 5 and Figure 6, respectively. For larger amplitude of harmonic film height the response in terms of both film pressure and total load are nonlinear with respect to the applied harmonic motion of the circular plate. Such nonlinearity is appropriately calculated by performing a transient analysis. Figure 6 also compares analytical and calculated time dependent values of the total load on the circular disc. The plot indicates that the calculated values agree well with the analytical values.
Figure 2: Film pressure vs radial distance from the center of the circular disc. The results computed by COMSOL are shown as the continuous curve and the analytical result is shown with green symbols.
Figure 3: Radial distribution of the film load on the wall
Figure 4: Film height variation vs time for different values of amplitude of harmonic film height.
Figure 5: Film pressure vs time for different values of amplitude of harmonic film height. For higher values of film height amplitude the film pressure varies nonlinearly with respect to the film height.
Figure 6: Total vertical load on the circular disc vs time. for different values of amplitude of harmonic film height. The results computed by COMSOL are shown by solid continuous curve and the analytical result is shown by dashed continuous curves. Similar to the film pressure, for higher values of film height amplitude the total load varies nonlinearly with respect to the film height.
Reference
1. B.J. Hamrock, S.R. Schmid, and B.O. Jacobson, Fundamentals of Fluid Film Lubrication, Marcel Dekker, 2004.
This model is based on the discussion entitled Parallel-Surface Bearing of infinite width in section 12.2 of the above reference.
Application Library path: MEMS_Module/Sensors/squeeze_film_disc
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Fluid Flow>Thin-Film Flow>Thin-Film Flow, Edge (tffs).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Frequency Domain.
6
Click  Done.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
r0
1[mm]
0.001 m
Disc radius
h0
10[um]
1E-5 m
Gap height
dh
dhND*h0
1E-7 m
Change in gap height
dhND
0.01
0.01
Fractional gap height change.
mu0
1e-5[Pa*s]
1E-5 Pa·s
Gas viscosity
f0
1000[Hz]
1000 Hz
Vibration frequency
Geometry 1
Line Segment 1 (ls1)
1
In the Geometry toolbar, click  More Primitives and choose Line Segment.
2
In the Settings window for Line Segment, locate the Starting Point section.
3
From the Specify list, choose Coordinates.
4
Locate the Endpoint section. From the Specify list, choose Coordinates.
5
In the r text field, type r0.
6
Click  Build All Objects.
7
Click the  Zoom Extents button in the Graphics toolbar.
Definitions
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 1 only.
Variables 1
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
In the table, enter the following settings:
Name
Expression
Unit
Description
Ftot
intop1(tffs.fwallz)
N
Total force on disc
Ftotan
-3*pi*mu0*vf*r0^4/(2*h0^3)
N
Analytic expression
vf
2*pi*f0*dh
m/s
Disc velocity
Ftotantime
-6*pi^2*mu0*f0*r0^4*dh*cos(2*pi*f0*t)/(2*(h0+dh*sin(2*pi*f0*t))^3)
N
Analytic expression
Analytic 1 (an1)
1
In the Definitions toolbar, click  Analytic.
2
In the Settings window for Analytic, type Pan in the Function name text field.
3
Locate the Definition section. In the Expression text field, type -6*pi*mu0*f0*dh*(r0^2-rf^2)/h0^3.
4
In the Arguments text field, type rf.
5
Locate the Units section. In the table, enter the following settings:
Argument
Unit
rf
m
6
In the Function text field, type Pa.
Materials
Material 1 (mat1)
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, locate the Material Contents section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Dynamic viscosity
mu
mu0
Pa·s
Basic
The density is not required because the modified Reynolds equation is solved.
Thin-Film Flow, Edge (tffs)
1
In the Model Builder window, under Component 1 (comp1) click Thin-Film Flow, Edge (tffs).
2
In the Settings window for Thin-Film Flow, Edge, locate the Physical Model section.
3
From the Equation type list, choose Modified Reynolds equation.
Fluid-Film Properties 1
Reverse the direction of the geometry normal. Note that this is necessary to ensure that the wall load acts in the vertical direction. See Figure 1 for more details on the orientation of the wall and base with respect to the geometric and reference surface normals. To view the orientation of the reference normal, show the default solver settings and then compute only to the dependent variables stage in the study sequence. Then plot the reference surface normal.
1
In the Model Builder window, under Component 1 (comp1)>Thin-Film Flow, Edge (tffs) click Fluid-Film Properties 1.
2
In the Settings window for Fluid-Film Properties, click to expand the Reference Surface Properties section.
3
From the Reference normal orientation list, choose Opposite direction to geometry normal.
4
Locate the Wall Properties section. In the hw1 text field, type h0.
5
From the uw list, choose None.
6
From the vw list, choose User defined. Specify the vector as
0
r
vf
z
Fluid-Film Properties 2
1
Right-click Component 1 (comp1)>Thin-Film Flow, Edge (tffs)>Fluid-Film Properties 1 and choose Duplicate.
2
Select Boundary 1 only.
3
In the Settings window for Fluid-Film Properties, locate the Wall Properties section.
4
From the uw list, choose User defined. Specify the associated vector as
0
r
dh*sin(2*pi*f0*t)
z
5
From the vw list, choose Calculate from wall displacement.
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Extra fine.
4
Click  Build All.
5
Click the  Zoom Extents button in the Graphics toolbar.
Study 1
Step 1: Frequency Domain
1
In the Model Builder window, under Study 1 click Step 1: Frequency Domain.
2
In the Settings window for Frequency Domain, locate the Study Settings section.
3
In the Frequencies text field, type 1000.
4
Locate the Physics and Variables Selection section. Select the Modify model configuration for study step check box.
5
In the tree, select Component 1 (Comp1)>Thin-Film Flow, Edge (Tffs)>Fluid-Film Properties 2.
6
Right-click and choose Disable.
7
In the tree, select Component 1 (Comp1).
8
In the Home toolbar, click  Compute.
Results
1D Plot Group 2
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, locate the Plot Settings section.
3
Select the x-axis label check box.
4
In the associated text field, type Radial distance (m).
5
Select the y-axis label check box.
6
In the associated text field, type Pressure (Pa).
7
Locate the Legend section. From the Position list, choose Upper left.
Line Graph 1
1
Right-click 1D Plot Group 2 and choose Line Graph.
2
Select Boundary 1 only.
3
In the Settings window for Line Graph, click to expand the Legends section.
4
Select the Show legends check box.
5
From the Legends list, choose Manual.
6
In the table, enter the following settings:
Legends
COMSOL
Line Graph 2
1
In the Model Builder window, right-click 1D Plot Group 2 and choose Line Graph.
2
Select Boundary 1 only.
3
In the Settings window for Line Graph, locate the y-Axis Data section.
4
In the Expression text field, type Pan(r).
5
Click to expand the Coloring and Style section. Find the Line style subsection. From the Line list, choose None.
6
Find the Line markers subsection. From the Marker list, choose Cycle.
7
Locate the Legends section. Select the Show legends check box.
8
From the Legends list, choose Manual.
9
In the table, enter the following settings:
Legends
Analytic Solution
10
In the 1D Plot Group 2 toolbar, click  Plot.
11
Click the  Zoom Extents button in the Graphics toolbar.
Radial Pressure
1
Right-click 1D Plot Group 2 and choose Rename.
2
In the Rename 1D Plot Group dialog box, type Radial Pressure in the New label text field.
3
Click OK.
2D Plot Group 3
In the Home toolbar, click  Add Plot Group and choose 2D Plot Group.
Arrow Line 1
1
Right-click 2D Plot Group 3 and choose Arrow Line.
2
In the Settings window for Arrow Line, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Thin-Film Flow, Edge>Fluid loads>tffs.fwallr,tffs.fwallz - Fluid load on wall.
3
Locate the Arrow Positioning section. In the Number of arrows text field, type 30.
4
Locate the Coloring and Style section. Select the Scale factor check box.
5
Use the slider to adjust the arrow length.
6
In the 2D Plot Group 3 toolbar, click  Plot.
7
Click the  Zoom Extents button in the Graphics toolbar.
Line 1
1
In the Model Builder window, right-click 2D Plot Group 3 and choose Line.
2
In the Settings window for Line, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Thin-Film Flow, Edge>Fluid loads>Fluid load on wall - N/m²>tffs.fwallz - Fluid load on wall, z component.
3
Locate the Coloring and Style section. From the Color table transformation list, choose Reverse.
4
In the 2D Plot Group 3 toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
Wall Load
1
Right-click 2D Plot Group 3 and choose Rename.
2
In the Rename 2D Plot Group dialog box, type Wall Load in the New label text field.
3
Click OK.
Global Evaluation 1
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Expressions section.
3
In the table, enter the following settings:
Expression
Unit
Description
abs(Ftot)
N
 
4
Click  Evaluate.
Global Evaluation 2
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, click Replace Expression in the upper-right corner of the Expressions section. From the menu, choose Component 1 (comp1)>Definitions>Variables>Ftotan - Analytic expression - N.
3
Clicknext to  Evaluate, then choose Table 1 - Global Evaluation 1.
Root
The agreement between the total force computed by COMSOL and the analytic expression is excellent.
Next add a study to solve the problem in the time domain.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies>Time Dependent.
4
Click Add Study in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Step 1: Time Dependent
1
In the Settings window for Time Dependent, locate the Study Settings section.
2
Click  Range.
3
In the Range dialog box, type 1/(40*f0) in the Step text field.
4
In the Stop text field, type 5/f0.
5
Click Replace.
6
In the Settings window for Time Dependent, click to expand the Study Extensions section.
7
Select the Auxiliary sweep check box.
8
Click  Add.
9
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
dhND (Fractional gap height change.)
0.01 0.1 0.2 0.3
 
10
In the Home toolbar, click  Compute.
Results
1D Plot Group 5
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
Point Graph 1
1
Right-click 1D Plot Group 5 and choose Point Graph.
2
Select Point 1 only.
3
In the Settings window for Point Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Thin-Film Flow, Edge>Wall and base properties>tffs.hw - Height of wall above reference plane - m.
4
In the 1D Plot Group 5 toolbar, click  Plot.
5
Click to expand the Legends section. Select the Show legends check box.
Wall height vs. time
1
In the Model Builder window, click 1D Plot Group 5.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
Select the Manual axis limits check box.
4
In the y minimum text field, type 0.
5
Locate the Legend section. From the Position list, choose Lower right.
6
In the 1D Plot Group 5 toolbar, click  Plot.
7
Right-click 1D Plot Group 5 and choose Rename.
8
In the Rename 1D Plot Group dialog box, type Wall height vs. time in the New label text field.
9
Click OK.
10
In the Wall height vs. time toolbar, click  Plot.
11
Click the  Zoom Extents button in the Graphics toolbar.
The wall height varies sinusoidally, as expected. Except in the case when dhND is 0.01, the height varies by a significant fraction of the gap. This means that the equation system is nonlinear.
1D Plot Group 6
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
Point Graph 1
1
Right-click 1D Plot Group 6 and choose Point Graph.
2
Select Point 1 only.
3
In the 1D Plot Group 6 toolbar, click  Plot.
Pressure vs. time
1
In the Model Builder window, right-click 1D Plot Group 6 and choose Rename.
2
In the Rename 1D Plot Group dialog box, type Pressure vs. time in the New label text field.
3
Click OK.
4
In the Pressure vs. time toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
Once the wall height varies by a significant fraction of the gap nonlinear effects become important. The pressure does not vary sinusoidally any more and it is no longer possible to model the equation straightforwardly in the frequency domain.
When the height variations are small direct comparison with the frequency domain solution is possible. Do this by comparing the maximum value of the pressure in the time domain with the amplitude of the pressure computed in the frequency domain.
Point Evaluation 1
1
In the Results toolbar, click  Point Evaluation.
2
Select Point 1 only.
3
In the Settings window for Point Evaluation, locate the Expressions section.
4
In the table, enter the following settings:
Expression
Unit
Description
abs(pfilm)
Pa
 
5
Click  Evaluate.
Point Evaluation 2
1
In the Results toolbar, click  Point Evaluation.
2
In the Settings window for Point Evaluation, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
4
From the Parameter selection (dhND) list, choose From list.
5
In the Parameter values (dhND) list, select 0.01.
6
Select Point 1 only.
7
Locate the Data Series Operation section. From the Transformation list, choose Minimum.
8
Click  Evaluate.
Table
1
Go to the Table window.
In this case the time and frequency domain analyses are in good agreement.
Results
1D Plot Group 7
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
Global 1
1
Right-click 1D Plot Group 7 and choose Global.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Definitions>Variables>Ftot - Total force on disc - N.
Global 2
1
In the Model Builder window, right-click 1D Plot Group 7 and choose Global.
2
In the Settings window for Global, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Definitions>Variables>Ftotantime - Analytic expression - N.
3
Click to expand the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dashed.
Total force on disc vs. time
1
In the Model Builder window, click 1D Plot Group 7.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
Select the Manual axis limits check box.
4
In the x minimum text field, type 0.003.
5
In the x maximum text field, type 0.005.
6
Locate the Legend section. From the Position list, choose Lower left.
7
Right-click 1D Plot Group 7 and choose Rename.
8
In the Rename 1D Plot Group dialog box, type Total force on disc vs. time in the New label text field.
9
Click OK.
10
In the Total force on disc vs. time toolbar, click  Plot.
11
Click the  Zoom Extents button in the Graphics toolbar.
The agreement between the COMSOL result and the analytic expression is good. The force acting on the disc remains periodic but has higher harmonic components at large displacements relative to the gap size.