COMSOL 6.4 - Squeeze Film Between Rectangular Porous and Nonporous Plates

Squeeze Film Between Rectangular Porous and Nonporous Plates

This example describes how to simulate the flow of a thin film of fluid in the gap between two rectangular plates, one of them with a porous facing, when the fluid is squeezed as a consequence of the relative motion between the plates. The model accounts for the ingress and egress of fluid between the porous material and the thin-film region. The simple geometry of the model allows for the comparison of the computed solution with analytic expressions in Ref. 1.

Model Definition

The model is set up with parameterized quantities in a unitless base system. It uses a thin-film flow approximation, described by the Reynolds equation, for the fluid motion in the gap between the plates. The fluid flow in the porous region is described by Darcy’s law. The coupling between the Thin-Film Flow interface and the Darcy’s Law interface is handled by the Thin-Film and Porous Media multiphysics coupling.

Geometry

The length of the rectangular plates, a, is taken to be of unit value. The width, b is controlled by the aspect ratio parameter, k = a/b. Using symmetry planes orthogonal to the x- and y-axes, only one quarter of the plate is modeled, as shown in Figure 1. The height of the porous facing, H, is suitably scaled with respect to the plate length using the parameter,

= H/a. A rigid body motion is considered for the nonporous plate and it is therefore not included in the model geometry. The relative motion between the porous and nonporous plates can be accounted for by the motion of the “wall” or “base” defined in the Thin-Film Flow interface.

Figure 1: Computational domain with symmetry boundaries (in yellow) and the thin-film region (in purple).

Physics Interface Settings

The Thin-Film Flow interface is set up to solve the Reynolds equation. The reference pressure, pref, is set to zero. The parameterized fluid density, ρ, and viscosity, η, are both set to unit value. The thin-film thickness is chosen to be, ht = 0.001, and set as the height of the wall above the reference plane. The squeeze action is induced by applying a wall velocity in the positive z direction of unit magnitude controlled by the parameter, vsq. Symmetry condition is applied to edges corresponding to the symmetry boundaries in the geometry. Border condition with zero pressure is applied to other edges as default.

The Darcy’s Law interface is set up with zero reference pressure and similar fluid properties to that in the Thin-Film Flow interface. The porosity, εp, is set to zero and the permeability,

, is chosen to be isotropic in nature, with a value controlled by the permeability parameter,

. The boundary conditions include a No Flow condition on the top face, a Symmetry condition on the symmetry boundaries, and a Pressure condition with zero prescribed value elsewhere.

The multiphysics coupling, added automatically, acts at the intersection of the two interfaces. It accounts for the fluid motion and the pressure balance in the region.

Meshing

A structured mesh with hexahedral elements is constructed using the Mapped and Swept mesh features. The porous domain thickness is meshed with four elements. The resulting mesh is shown in Figure 2.

Figure 2: Computational domain with structured, hexahedral mesh.

Results and Discussion

The analysis presented in Ref. 1 provides the expression for the dimensionless pressure at the intersection of the physics interfaces. It is given in the form of infinite series as

(1)

where

(2)

(3)

(4)

(5)

The dimensionless load carrying capacity of the squeeze film is given as

(6)

In the comparison plots presented in the following paragraphs, the values of m and n are taken to be equal to 29, as suggested in Ref. 1. This ensures accurate pressure and load capacity values up to five significant digits.

Figure 3 shows pressure distribution for different values of ψ on the intersecting boundary of the physics interfaces at the location y = b/2, for values k = 1 and

= 0.02. The maximum pressure rapidly decreases with an increase in the permeability parameter. For values of ψ < 0.001, the behavior is close to the nonporous limit as reported in Ref. 1.

Figure 3: Pressure distribution for various values of ψ.

Figure 4 shows the load capacity for different values of k at the intersecting boundary of the physics interfaces for the value of

= 0.02. The load-carrying capacity generally decreases with an increase in the aspect ratio, as such plates cannot contain the pressure generated due to increased side leakage. In fact, square plates provide the highest load-bearing capacity for a given quadrilateral area of coverage.

Figure 4: Load capacity for various values of k.

Figure 5 shows the loss of load capacity with increasing aspect ratio and permeability parameter.

Figure 5: Load capacity for various values of ψ.

Figure 6 shows load capacity for different values of

at the intersection of the physics interfaces for the values of k = 1 and ht = 0.001. The load capacity is large and shows little variation with changes in the porous domain thickness for low permeabilities.With an increase in the permeability and porous domain thickness, however, the load capacity rapidly deteriorates.

Figure 6: Load capacity for various values of

All comparison plots show excellent agreement between computed and analytic solutions.

Figure 7 shows the pressure developed at the intersection of the physics interfaces and the subsequent percolation flow in the porous domain as a consequence of the squeezing action of the plates on the fluid in the thin-film region for the case k = 1,

= 0.02, and ψ = 0.01.

Figure 7: Darcy pressure and velocity streamlines.

Reference

1. H. Wu, “An Analysis of the Squeeze Film Between Porous Rectangular Plates,” J. Lubr. Technol., vol. 94, no. 1, pp. 64–68, 1972.

CFD_Module/Verification_Examples/squeeze_film_porous

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Choose a unitless base system.

Root

In the window, click the root node.

In the root node’s window, locate the section.

From the list, choose .

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file squeeze_film_porous_parameters.txt.

Define analytic expressions of the solution. The expressions are available in Equation 1–Equation 6.

Definitions

Analytic 1 (an1)

In the toolbar, click

In the window for , type a_m in the text field.

Locate the section. In the text field, type m*pi/a.

In the text field, type m, a.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	m	0	29	0

Analytic 2 (an2)

In the toolbar, click

In the window for , type b_n in the text field.

Locate the section. In the text field, type n*pi/b.

In the text field, type n, b.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	n	0	29	0

Analytic 3 (an3)

In the toolbar, click

In the window for , type gamma_mn in the text field.

Locate the section. In the text field, type pi*sqrt(m^2+(k*n)^2)/a.

In the text field, type m, n, k, a.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	m	0	29	0
√	n	0	29	0
√	k	0	8	0

Analytic 4 (an4)

In the toolbar, click

In the window for , type G_mn in the text field.

Locate the section. In the text field, type 12*(exp(2*gamma_mn*H)-1)/(H*(exp(2*gamma_mn*H)+1)/a).

In the text field, type gamma_mn, H, a.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	gamma_mn	0	1e6	0

Analytic 5 (an5)

In the toolbar, click

In the window for , type p_summand in the text field.

Locate the section. In the text field, type (sin(a_m((2*m-1),a)*x)*sin(b_n((2*n-1),b)*y))/((2*m-1)*(2*n-1)*sqrt((2*m-1)^2+(k*(2*n-1))^2)*(pi*sqrt((2*m-1)^2+(k*(2*n-1))^2)+(psi*G_mn(gamma_mn((2*m-1),(2*n-1),k,a),H,a)))).

In the text field, type m, n, x, y, a, b, psi, k, H.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	m	0	29	0
√	n	0	29	0
√	x	0	a	0

Analytic 6 (an6)

In the toolbar, click

In the window for , type W_summand_psi in the text field.

Locate the section. In the text field, type 1/(((2*m-1)*(2*n-1))^2*sqrt((2*m-1)^2+(k*(2*n-1))^2)*((pi*sqrt((2*m-1)^2+(k*(2*n-1))^2))+(psi*G_mn(gamma_mn((2*m-1),(2*n-1),k,a),H,a)))).

In the text field, type m, n, a, psi, k, H.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	m	0	29	0
√	n	0	29	0

Analytic 7 (an7)

In the toolbar, click

In the window for , type W_summand_phi in the text field.

Locate the section. In the text field, type 1/(((2*m-1)*(2*n-1))^2*sqrt((2*m-1)^2+(k*(2*n-1))^2)*((pi*sqrt((2*m-1)^2+(k*(2*n-1))^2))+(phi*H/ht^3*G_mn(gamma_mn((2*m-1),(2*n-1),k,a),H,a)))).

In the text field, type m, n, a, phi, k, H.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	m	0	29	0
√	n	0	29	0

Geometry 1

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type a/2.

In the text field, type b/2.

In the text field, type H.

Form Union (fin)

In the toolbar, click

Create selections for the thin-film domain and symmetry boundaries.

Definitions

Thin-Film Boundary

In the toolbar, click

In the window for , type Thin-Film Boundary in the text field.

Locate the section. From the list, choose .

Click

In the dialog, type 3 in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Symmetry Boundaries

In the toolbar, click

In the window for , type Symmetry Boundaries in the text field.

Locate the section. From the list, choose .

Click

In the dialog, type 5-6 in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Click the

button in the toolbar.

Click the

button in the toolbar.

Next, specify parameters for the Thin-Film Flow interface: Set the reference pressure to zero and reverse the reference normal to point in the positive z direction.

Thin-Film Flow (tff)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the pref text field, type 0.

Fluid-Film Properties 1

In the window, under > click .

In the window for , locate the section.

From the T list, choose .

Click to expand the section. From the list, choose .

Locate the section. In the hw1 text field, type ht.

From the vw list, choose . Specify the vector as


sqvel	z

Locate the section. From the μ list, choose . In the associated text field, type eta.

From the ρ list, choose . In the associated text field, type rho.

Symmetry 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog, type 7, 10 in the text field.

Click .

Now specify parameters for the Darcy’s Law interface, including a zero reference pressure.

Darcy’s Law (dl)

In the window, under click .

In the window for , locate the section.

In the pref text field, type 0.

Fluid 1

In the window, under > > click .

In the window for , locate the section.

From the ρ list, choose . In the associated text field, type rho.

From the μ list, choose . In the associated text field, type eta.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the εp list, choose . From the κ list, choose . In the associated text field, type phi.

Symmetry 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Pressure 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog, type 1, 2 in the text field.

Click .

Create a structured hexahedral mesh.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Mapped 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Swept 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 4.

Free Tetrahedral 1

In the window, under > right-click and choose .

Right-click and choose .

Click the

button in the toolbar.

Study 1

In the window, click .

In the window for , locate the section.

Clear the checkbox.

Add a Parametric Sweep step for computing the Stationary solution for specified parameter combinations.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
k (Ratio of plate length to width)	81^range(1,5) 41^range(1,5) 2*1^range(1,5) 1^range(1,5) 1^range(1,15)

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
psi (Permeability parameter)	10^{range(0,-1,-4)} 10^{range(0,-1,-4)} 10^{range(0,-1,-4)} 10^{range(0,-1,-4)} range(0.1,0.1,0.5) range(0.01,0.01,0.05) range(0.001,0.001,0.005)

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
Hbar (Ratio of porous domain thickness to plate length)	0.02*1^range(1,20) range(0.02,0.02,0.1) range(0.02,0.02,0.1) range(0.02,0.02,0.1)