Forchheimer Flow

This is a tutorial of the coupling between flow of a fluid in an open channel and a porous block attached to one of the channel walls. The flow is described by the Navier-Stokes equation in the free region and a Forchheimer-corrected version of the Brinkman equations in the porous region.

Figure 1: Depiction of the modeling geometry and domain. The 3D geometry can be reduced to a 2D representation assuming that changes through the thickness are negligible.

The coupling of free media flow with porous media flow is common in the fields of earth science and chemical engineering. Perhaps the most frequent way to deal with coupled free and porous media flow is to incorporate Darcy’s law adjacent to Navier-Stokes because it is usually numerically easy to solve. However, this approach does not account for viscous effects arising from the free media flow, which may still be important in the region close to the free-porous structure interface. Depending on the pore size and pore distribution, but also the fluid’s properties, it can therefore be an oversimplification to employ Darcy’s law. The Brinkman equations account for momentum transport by macroscopic viscous effects as well as pressure gradients (stemming from microscopic shear effects inside each pore channel) and can be considered an extension of Darcy’s Law.

Still, the Brinkman equations assume laminar flow. Looking at processes in relatively open structures, like gas flow through packed beds, there is also a turbulent contribution to the resistance to flow. In those cases, an additional term accounts for the turbulent contribution to the resistance to flow in the porous domain. The Forchheimer equation (also accredited to Ergun) is widely used to predict pressure drops in packed beds. This equation can generally be written as

The left-hand side is the pressure drop per unit length of traveled distance through the bed. The first term on the right-hand side represents the Blake-Kozeny equation for laminar flow. The pressure drop depends linearly on the average linear velocity u for laminar flow, corresponding to Darcy flow. The second term is from the purely turbulent Burke-Plummer equation where the pressure drop is proportional to the square of the velocity. Description of an intermediate flow, where both the laminar and turbulent effects are important, requires the two-term Forchheimer equation. The coefficients α1 and α2 are functions of porosity, viscosity, average pore diameter, and fluid density.

Model Definition

Figure 2 below shows the example domain and notations for the boundary conditions.

Figure 2: Modeled domain and boundary notations. Flow enters at the bottom and leaves at the top. The region of porous structure is not as long as the free channel.

Flow in the free channel is described by the stationary, incompressible Navier-Stokes equations:

(1)

where μ denotes the dynamic viscosity (Pa·s), u refers to the velocity in the open channel (m/s), ρ is the fluid’s density (kg/m3), and p is the pressure (Pa). In the porous domain, the Brinkman equations with Forchheimer correction describe the flow:

(2)

Here k denotes the permeability of the porous medium (m2), εp is the porosity (dimensionless), and the dimensionless friction coefficient is (Ref. 1)

As Equation 1 and Equation 2 reveal, the momentum transport equations are closely related. The term on the left-hand side of the Navier-Stokes formulation corresponds to momentum transferred by convection in free flow. The Brinkman formulation replaces this term by a contribution associated with the drag force experienced by the fluid flowing through a porous medium. In addition, the last term in the right-hand side of Equation 2 presents the Forchheimer correction for turbulent drag contributions.

In COMSOL Multiphysics, it is easy to set up and solve such a coupled flow regime. The implementation of the extra drag is done with a Forchheimer coefficient (kg/m4) equal to

The boundary conditions are

where the pressure level at the outlet is used as a reference value.

The following table lists the input data for the example:


Property	Value	Description
μ	10-3 kg/(m·s)	Dynamic viscosity
ρ	1000 kg/m3	Density
k	10-7 m2	Permeability
εp	0.4	Porosity
v0	2 cm/s	Inlet velocity

Results and Discussion

Figure 3 shows the velocity field in the open channel and porous structure. The plot reveals that there are slight disturbances in the velocity at the porous wall, which suggests momentum transport by viscous effects.

Figure 3: Velocity field without Forchheimer correction.

Figure 4 shows the corresponding field in the porous medium, as modeled by the Forchheimer-corrected Brinkman equation.

Figure 4: Velocity field in the both free flow and porous medium.

Figure 5 contains a cross-sectional velocity plot. It shows that without the Forchheimer correction, the resistance to flow is underestimated in the porous domain. The added correction gives a solution with slower flow in the porous domain and a faster flow in the free domain, due to incompressibility.

Figure 5: Cross section of the velocity field (velocity magnitude) in the middle of the modeling domain with and without the Forchheimer correction (stars and circles markers, respectively).

Further postprocessing would show that the shear rate perpendicular to the flow is also continuous. This implies that there is significant viscous momentum transfer across the interface and into the porous material, a transport that is not accounted for by Darcy’s law.

Notes About the COMSOL Implementation

To implement the Forchheimer pressure drop relation in a differential equation framework, this example uses an approach suggested in Ref. 1 in which the Brinkman momentum balance is amended with a Forchheimer term. The system studied in this example is that of a 2D cross section of a rectangular channel with a porous layer attached to one of its walls. Flow enters the volume with a uniform velocity profile and develops throughout the length of the channel.

Reference

1. A. Amiri and K. Vafai, “Transient Analysis of Incompressible Flow Through a Packed Bed,” Int. J. Heat and Mass Transfer, vol. 41, pp. 4259–4279, 1998.

Subsurface_Flow_Module/Fluid_Flow/forchheimer_flow

Modeling Instructions

From the menu, choose .

New

In the window, click .

Model Wizard

In the window, click .

In the tree, select .

Click .

In the tree, select .

Click .

Geometry 1

Rectangle 1 (r1)

On the toolbar, click and choose .

In the window for , locate the section.

In the text field, type 1e-3.

In the text field, type 6e-3.

Locate the section. In the text field, type -3e-3.

Rectangle 2 (r2)

On the toolbar, click and choose .

In the window for , locate the section.

In the text field, type 1e-3.

In the text field, type 8e-3.

Locate the section. In the text field, type -1e-3.

In the text field, type -4e-3.

On the toolbar, click .

Global Definitions

Parameters

On the toolbar, click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
v0	2[cm/s]	0.02 m/s	Inlet velocity
eps_p	0.4	0.4	Porosity
Cf	1.75/sqrt(150*eps_p^3)	0.56481	Friction coefficient
fs	1	1	Switch for Forchheimer terms

Materials

Add blank materials for the fluid and porous matrix. You will specify their properties after the set up of the physics interface.

Material 1 (mat1)

In the window, under right-click and choose .

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

Material 2 (mat2)

Right-click and choose .

In the window for , type Porous Matrix in the text field.

Laminar Flow (spf)

In the window, under click .

In the window for , locate the section.

Select the check box.

Fluid and Matrix Properties 1

On the toolbar, click and choose .

Select Domain 2 only.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Now, COMSOL recognizes which material properties are required to solve this model.

Materials

Fluid (mat1)

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Name	Value	Unit	Property group
Density	rho	1000[kg/m^3]	kg/m³	Basic
Dynamic viscosity	mu	1e-3[Pa*s]	Pa·s	Basic

Porous Matrix (mat2)

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Name	Value	Unit	Property group
Porosity	epsilon	eps_p	1	Basic
Permeability	kappa	1e-7[m^2]	m²	Basic

Laminar Flow (spf)

Fluid and Matrix Properties 1

In the window, under click .

Forchheimer Drag 1

On the toolbar, click and choose .

In the window for , locate the section.

In the βF text field, type fs*Cf*eps_p*spf.rho/sqrt(spf.kappaxx).

Recall that fs is a parameter switching on and off the Forchheimer drag. Moreover, spf.rho is the density variable for the Laminar Flow interface and spf.kappabr the permeability. To find these names, click in the Model Builder window’s toolbar and choose . Then go to the nodes for the Fluid and Matrix Properties.

Inlet 1

On the toolbar, click and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the U0 text field, type v0.

Outlet 1

On the toolbar, click and choose .

Select Boundary 3 only.

Mesh 1

Size