PDF

Forchheimer Flow
Introduction
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:
μ
10-3 kg/(m·s)
ρ
10-7 m2
v0
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.
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.
Application Library path: Subsurface_Flow_Module/Fluid_Flow/forchheimer_flow
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.
2
In the Select Physics tree, select Fluid Flow>Porous Media and Subsurface Flow>Free and Porous Media Flow (fp).
3
Click  Add Physics.
4
Click  Study.
5
In the Select Study tree, select General Studies>Stationary.
6
Geometry 1
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape 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
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape 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.
6
In the y text field, type -4e-3.
7
In the Geometry toolbar, click  Build All.
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
Materials
Add blank materials for the fluid and porous matrix. You will specify their properties after the setup of the physics interface.
Fluid
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.
Porous Matrix
1
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 and Matrix Properties 1
1
In the Model Builder window, under Component 1 (comp1) right-click Free and Porous Media Flow (fp) and choose Fluid and Matrix Properties.
2
3
In the Settings window for Fluid and Matrix Properties, locate the Fluid Properties section.
4
From the Fluid material list, choose Fluid (mat1).
5
Locate the Porous Matrix Properties section. From the Porous material list, choose Porous Matrix (mat2).
Now, COMSOL Multiphysics recognizes which material properties are required to solve this model.
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
Porous Matrix (mat2)
1
In the Model Builder window, click Porous Matrix (mat2).
2
In the Settings window for Material, locate the Material Contents section.
3
Free and Porous Media Flow (fp)
Fluid and Matrix Properties 1
1
In the Model Builder window, under Component 1 (comp1)>Free and Porous Media Flow (fp) click Fluid and Matrix Properties 1.
2
In the Settings window for Fluid and Matrix Properties, locate the Porous Matrix Properties section.
3
From the Permeability model list, choose Non-Darcian.
4
In the cF text field, type fs*Cf.
Recall that fs is a parameter switching on and off the Forchheimer drag.
Inlet 1
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
3
In the Settings window for Inlet, locate the Velocity section.
4
In the U0 text field, type v0.
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Mesh 1
Free Triangular 1
In the Mesh toolbar, click  Free Triangular.
Size
1
In the Model Builder window, click Size.
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  Build All.
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
Select the Auxiliary sweep check box.
4
5
6
In the Home toolbar, click  Compute.
Results
Velocity (fp)
Visualize the velocity fields as in Figure 3 and Figure 4.
Arrow Surface 1
1
Right-click Velocity (fp) and choose Arrow Surface.
2
In the Velocity (fp) toolbar, click  Plot.
3
Click the  Zoom Extents button in the Graphics toolbar.
Velocity (fp)
1
In the Model Builder window, click Velocity (fp).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Parameter value (fs) list, choose 0.
4
In the Velocity (fp) toolbar, click  Plot.
Create the plot from Figure 5.
Cut Line 2D 1
1
In 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.
Velocity magnitude
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Velocity magnitude in the Label text field.
Line Graph 1
1
Right-click Velocity magnitude and choose Line Graph.
2
In the Settings window for Line Graph, locate the Data section.
3
From the Dataset list, choose Cut Line 2D 1.
4
In the Velocity magnitude toolbar, click  Plot.
5
Click to expand the Coloring and Style section. Find the Line markers subsection. From the Marker list, choose Cycle.
6
Click to expand the Legends section. Select the Show legends check box.
7
From the Legends list, choose Manual.
8
9
In the Velocity magnitude toolbar, click  Plot.