Solved with COMSOL Multiphysics 5.

Forchheimer Flow
Introduction
This is a tutorial model 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.
Free flow Porous structure

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

1 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

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

Δp
------- = α 1 u + α 2 u 2
L

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 model domain and notations for the boundary conditions.

Outlet

No slip No slip

Navier-Stokes
Brinkman boundary

Inlet

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:

2 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

T
ρ ( u ⋅ ∇ )u = – ∇ ⋅ [ – p I + μ ( ∇u + ( ∇u ) ) ] (1)
∇⋅u = 0

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:

μ μ- T ρε p C f
--- u = ∇ ⋅ – p I + ---- ( ∇u + ( ∇u ) ) – --------------- u u
k εp k (2)
∇⋅u = 0

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

1,75
C f = --------------------
3
150ε p

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, the predefined multiphysics interface Free and Porous

Media Flow makes it 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

ρε p C f
β F = ---------------
k

The boundary conditions are

u = ( u 0, v 0 ) Inlet: velocity
u = 0 Wall:no slip
T Outlet: Pressure, no viscous stress
p = 0, μ ( ∇u + ( ∇u ) ) = 0

3 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

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

The following table lists the input data for the model:

PROPERTY VALUE DESCRIPTION

-3
μ 10 kg/(m·s) Dynamic viscosity
ρ 1000 kg/m3 Density
-7 2
k 10 m 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 in the free channel.

4 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

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.

The velocity profile along a horizontal cross section in the middle of the channel is
shown in Figure 5. The plot confirms the existence of velocity gradients in the porous
structure in that the velocity is continuous across the interface between the free
channel and the porous structure. 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. Coupling Darcy’s law to the
Navier-Stokes equations in this example, would give a substantially different profile in
which the velocity in the free channel falls to a much smaller value in the free channel.
Furthermore, a constant velocity equal to the flat part of the curve in Figure 5 would
describe all transport in the porous media.

Figure 5 contains a cross-sectional velocity plot. It shows that without the

Forchheimer correction, the resistance to flow is underestimated in the porous

5 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

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 (squares and triangle
markers, respectively).

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.

6 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

Application Library path: Subsurface_Flow_Module/Fluid_Flow/

forchheimer_flow

Modeling Instructions
From the File menu, choose New.

NEW
1 In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select physics tree, select Fluid Flow>Porous Media and Subsurface Flow>Free
and Porous Media Flow (fp).
3 Click Add.
4 Click Study.
5 In the Select study tree, select Preset Studies>Stationary.
6 Click Done.

GEOMETRY 1

Rectangle 1 (r1)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size section.
3 In the Width text field, type 1e-3.
4 In the Height text field, type 6e-3.
5 Locate the Position section. In the y text field, type -3e-3.

Rectangle 2 (r2)
1 On the Geometry toolbar, click Primitives and choose Rectangle.
2 In the Settings window for Rectangle, locate the Size section.
3 In the Width text field, type 1e-3.
4 In the Height text field, type 8e-3.
5 Locate the Position section. In the x text field, type -1e-3.

7 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

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

7 Click the Zoom Extents button on the Graphics toolbar.

GLOBAL DEFINITIONS

Parameters
1 On the Home toolbar, click Parameters.
2 In the Settings window for Parameters, locate the Parameters section.
3 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.5648 Friction coefficient
fs 1 1 Switch Forchheimer term

MATERIALS
Add materials for the fluid and the porous matrix. You will specify their properties
shortly.

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, type Fluid in the Label text field.

Material 2 (mat2)
1 In the Model Builder window, right-click Materials and choose Blank Material.
2 In the Settings window for Material, type Porous Matrix in the Label text field.

FREE AND POROUS MEDIA FLOW (FP)

Fluid Properties 1
1 In the Model Builder window, expand the Free and Porous Media Flow (fp) node, then
click Fluid Properties 1.
2 In the Settings window for Fluid Properties, locate the Fluid Properties section.
3 From the Fluid material list, choose Fluid (mat1).

Porous Matrix Properties 1

1 On the Physics toolbar, click Domains and choose Porous Matrix Properties.

8 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

2 Select Domain 2 only.

3 In the Settings window for Porous Matrix Properties, locate the Porous Matrix
Properties section.
4 From the Porous material list, choose Porous Matrix (mat2).
Now go back and specify the material properties.

MATERIALS

Fluid (mat1)
1 In the Model Builder window, under Component 1 (comp1)>Materials click Fluid
(mat1).
2 In the Settings window for Material, locate the Material Contents section.
3 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)

1 In the Model Builder window, under Component 1 (comp1)>Materials click Porous
Matrix (mat2).
2 In the Settings window for Material, locate the Material Contents section.
3 In the table, enter the following settings:

Property Name Value Unit Property group

Permeability kappa 1e-7[m^2] m² Basic
Porosity epsilon eps_p 1 Basic

FREE A N D PO ROUS MEDIA FLOW (FP)

Forchheimer Drag 1
1 On the Physics toolbar, click Attributes and choose Forchheimer Drag.
2 In the Settings window for Forchheimer Drag, locate the Forchheimer Drag section.
3 In the βF text field, type fs*Cf*eps_p*fp.rho/sqrt(fp.kappabr).
Recall that fs is a parameter switching on and off the Forchheimer drag. Moreover,
fp.rho is the density variable for the Free and Porous Media Flow interface and
fp.kappabr the permeability. To find these names, click Show in the Model Builder

9 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

window's toolbar and choose Equation View, and then go to the Equation View
nodes for the Fluid Properties and Porous Matrix Properties nodes, respectively.

Inlet 1
1 On the Physics toolbar, click Boundaries and choose Inlet.
2 Select Boundary 2 only.
3 In the Settings window for Inlet, locate the Velocity section.
4 In the U0 text field, type v0.

Outlet 1
1 On the Physics toolbar, click Boundaries and choose Outlet.
2 Select Boundary 3 only.

MESH 1

Size
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and
choose Free Triangular.
2 In the Settings window for Size, locate the Element Size section.
3 Click the Custom button.
4 Locate the Element Size Parameters section. In the Maximum element size text field,
type 1e-4.
5 Click the Build All button.

STUDY 1

Step 1: Stationary
1 In the Model Builder window, under Study 1 click Step 1: Stationary.
2 In the Settings window for Stationary, click to expand the Study extensions section.
3 Locate the Study Extensions section. Select the Auxiliary sweep check box.
4 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

fs 0 1

6 On the Home toolbar, click Compute.

10 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

RESULTS

Velocity (fp)
1 In the Model Builder window, under Results right-click Velocity (fp) and choose Arrow
Surface.
2 On the Velocity (fp) toolbar, click Plot.
3 Click the Zoom Extents button on the Graphics toolbar.
4 In the Model Builder window, click Velocity (fp).
5 In the Settings window for 2D Plot Group, locate the Data section.
6 From the Parameter value (fs) list, choose 0.
7 On the Velocity (fp) toolbar, click Plot.

Data Sets
1 On the Results toolbar, click Cut Line 2D.
2 In the Settings window for Cut Line 2D, locate the Line Data section.
3 In row Point 1, set x to -1e3.
4 In row Point 2, set x to 1e3.

1D Plot Group 3
1 On the Results toolbar, click 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Velocity magnitude in the Label
text field.

Velocity magnitude
1 On the Velocity magnitude toolbar, click Line Graph.
2 In the Settings window for Line Graph, locate the Data section.
3 From the Data set list, choose Cut Line 2D 1.
4 On the Velocity magnitude toolbar, click Plot.
5 Click to expand the Coloring and style section. Locate the Coloring and Style section.
Find the Line markers subsection. From the Marker list, choose Cycle.
6 On the Velocity magnitude toolbar, click Plot.

11 | FORCHHEIMER FLOW
 Solved with COMSOL Multiphysics 5.1

12 | FORCHHEIMER FLOW