Created in COMSOL Multiphysics 5.

Free Convection in a Porous Medium

This model is licensed under the COMSOL Software License Agreement 5.5.
All trademarks are the property of their respective owners. See www.comsol.com/trademarks.
Introduction
This example describes subsurface flow in porous media driven by density variations that
result from temperature changes. The example comes from Hossain and Wilson (Ref. 1),
who use a specialized in-house code to solve this free-convection problem. This COMSOL
Multiphysics example reproduces their work using the Brinkman Equations interface and
the Heat Transfer in Porous Media interface. The results of this model match those of the
published study.

Model Definition
The following figure gives the example geometry. Water in a porous medium layer can
move within the layer but not exit from it. Temperatures vary from high to low along the
outer edges. Initially the water is stagnant, but temperature gradients alter the fluid density
to the degree that buoyant flow occurs. The problem statement specifies that the flow is
steady state.
Tc

Th - (Th - Tc)s
Tc

Tc - (Tc - Th)s
Th

Figure 1: Domain geometry and boundary conditions for the heat balance in the free-
convection problem. Th is a higher temperature than Tc, while s is a variable that represents
the relative length of a boundary segment and goes from 0 to 1 along the segment.

Model this free-convection problem by introducing a Boussinesq buoyancy term to

Brinkman’s momentum equation, and then linking the resulting fluid velocities to the
Heat Transfer in Porous Media interface.

The Boussinesq buoyancy term that appears on the right-hand side of the momentum
equation accounts for the lifting force due to thermal expansion

2 | FREE CONVECTION IN A POROUS MEDIUM

 
--- u + p –   
---  u +  u  T  = g p  T – T c  (1)
 

u = 0.

In these expressions, T represents temperature, while Tc is a reference temperature, g

denotes the gravity acceleration,  gives the fluid density at the reference temperature,  is
the porosity, and p is the fluid’s coefficient of volumetric thermal expansion.

The heat balance comes from the heat transfer equation

C L u  T –    k eq T  = 0 (2)

where keq denotes the effective thermal conductivity of the fluid-solid mixture, and CL is
the fluid’s heat capacity at constant pressure.

The boundary conditions for the Brinkman equations are all no-slip conditions. Using
only velocity boundaries gives no information on the pressure within the domain, which
means that the example produces estimates of the pressure change instead of the pressure
field. However, without any seed information on pressure, the problem is unlikely to
converge. The remedy is to arbitrarily fix the pressure at a point in the example using a
point constraint. The boundary conditions for the Heat Transfer interface are the series of
fixed temperatures shown in Figure 1.

Implementation: Initial Conditions for Boussinesq Approximation

The simple statements in Equation 1 and Equation 2 produce a strong nonlinear problem
that represents a difficult convergence task for most nonlinear solvers. To ease the
numerical difficulties, let the coefficient of volumetric thermal expansion p increase
gradually, raising the Rayleigh number of the experiment. When p  0, the momentum
and temperature equations are uncoupled, so the example converges easily. Then increase
p, using the previous solution as the initial guess for the next parametric step, and so on,
until reaching a Rayleigh number of 105. The iteration protocol is an easy process with the
parametric solver in COMSOL Multiphysics.

Results
This example reproduces a model reported by Hossain and Wilson (Ref. 1). After
extracting the input data from the paper, the author constructed the example in less than

3 | FREE CONVECTION IN A POROUS MEDIUM

an hour, including all the steps from geometry input to postprocessing of the results.
Figure 2 shows the temperature distribution throughout the porous slice.

Figure 2: Temperature in a porous structure subjected to temperature gradients and

subsequent free convection. The COMSOL Multiphysics simulation is in excellent agreement
with published results from Ref. 1.

Figure 3 gives the COMSOL Multiphysics solution for the flow field.

4 | FREE CONVECTION IN A POROUS MEDIUM

Figure 3: Velocity field for a Rayleigh number of Ra = 105.

Reference
1. M. Anwar Hossain and M. Wilson, “Natural Convection Flow in a Fluid-saturated
Porous Medium Enclosed by Non-isothermal Walls with Heat Generation,” Int. J.
Therm. Sci., vol. 41, pp. 447–454, 2002.

Application Library path: Porous_Media_Flow_Module/Heat_Transfer/

convection_porous_medium

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

5 | FREE CONVECTION IN A POROUS MEDIUM

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select Physics tree, select Fluid Flow>Porous Media and Subsurface Flow>
Brinkman Equations (br).
3 Click Add.
4 In the Select Physics tree, select Heat Transfer>Heat Transfer in Porous Media (ht).
5 Click Add.
6 Click Study.
7 In the Select Study tree, select General Studies>Stationary.
8 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

rho 1000[kg/m^3] 1000 kg/m³ Fluid density
mu 0.001[Pa*s] 0.001 Pa·s Fluid viscosity
alphap 1e-6[1/K] 1E-6 1/K Fluid volumetric
thermal expansion
kth 6[W/(m*K)] 6 W/(m·K) Fluid thermal
conductivity
gamma 1 1 Fluid ratio of
specific heat
Cp 4200[J/(kg*K)] 4200 J/(kg·K) Fluid heat capacity
at constant pressure
epsilon 0.4 0.4 Porosity
kappa 1e-3[m^2] 0.001 m² Permeability
p0 1[atm] 1.0133E5 Pa Reference pressure
Tc 20[degC] 293.15 K Reference
temperature
Th 42[degC] 315.15 K High temperature
L 0.1[m] 0.1 m Length scale

6 | FREE CONVECTION IN A POROUS MEDIUM

Name Expression Value Description
Pr mu*Cp/kth 0.7 Prandtl number
Ra Cp*rho^2*g_const* 1.5102E5 Rayleigh number
alphap*(Th-Tc)*
L^3/(kth*mu)

GEOMETRY 1

Square 1 (sq1)
1 In the Geometry toolbar, click Square.
2 In the Settings window for Square, locate the Size section.
3 In the Side length text field, type L.
4 Click Build Selected.

Point 1 (pt1)
1 In the Geometry toolbar, click Point.
2 In the Settings window for Point, locate the Point section.
3 In the x text field, type L.
4 In the y text field, type L/10.
5 Click Build All Objects.

BRINKMAN EQUATIONS (BR)

Fluid and Matrix Properties 1

1 In the Model Builder window, under Component 1 (comp1)>Brinkman Equations (br) click
Fluid and Matrix Properties 1.
2 In the Settings window for Fluid and Matrix Properties, locate the Fluid Properties section.
3 From the  list, choose User defined. In the associated text field, type rho.
4 From the  list, choose User defined. In the associated text field, type mu.
5 Locate the Porous Matrix Properties section. From the p list, choose User defined. In the
associated text field, type epsilon.
6 From the  list, choose User defined. In the associated text field, type kappa.

Volume Force 1
1 In the Physics toolbar, click Domains and choose Volume Force.
Set up the Boussinesq buoyance term according to Equation 1.
2 In the Settings window for Volume Force, locate the Domain Selection section.

7 | FREE CONVECTION IN A POROUS MEDIUM

3 From the Selection list, choose All domains.
4 Locate the Volume Force section. Specify the F vector as

0 x
rho*g_const*alphap*(T-Tc) y

Pressure Point Constraint 1

1 In the Physics toolbar, click Points and choose Pressure Point Constraint.
2 Select Point 4 only.

HEAT TRANSFER IN POROUS MEDIA (HT)

Use quadratic elements for the discretization of the temperature field to improve accuracy
for this strongly coupled problem.

1 In the Model Builder window, under Component 1 (comp1) click

Heat Transfer in Porous Media (ht).
2 In the Settings window for Heat Transfer in Porous Media, click to expand the
Discretization section.
3 From the Temperature list, choose Quadratic Lagrange.
4 Click to collapse the Discretization section.

Porous Medium 1
1 In the Model Builder window, under Component 1 (comp1)>
Heat Transfer in Porous Media (ht) click Porous Medium 1.
2 In the Settings window for Porous Medium, locate the Model Input section.
3 From the pA list, choose User defined. In the associated text field, type p0.
4 Locate the Heat Convection section. From the u list, choose Velocity field (br).
5 Locate the Heat Conduction, Fluid section. From the k list, choose User defined. In the
associated text field, type kth.
6 Locate the Thermodynamics, Fluid section. From the  list, choose User defined. In the
associated text field, type rho.
7 From the Cp list, choose User defined. In the associated text field, type Cp.
8 From the  list, choose User defined. In the associated text field, type gamma.
9 Locate the Immobile Solids section. In the p text field, type 1-epsilon.
10 Locate the Heat Conduction, Porous Matrix section. From the kp list, choose User defined.
Locate the Thermodynamics, Porous Matrix section. From the p list, choose User defined.
From the Cp, p list, choose User defined.

8 | FREE CONVECTION IN A POROUS MEDIUM

Initial Values 1
1 In the Model Builder window, click Initial Values 1.
2 In the Settings window for Initial Values, locate the Initial Values section.
3 In the T text field, type Tc.

Temperature 1
1 In the Physics toolbar, click Boundaries and choose Temperature.
2 Select Boundary 2 only.
3 In the Settings window for Temperature, locate the Temperature section.
4 In the T0 text field, type Th.

Temperature 2
1 In the Physics toolbar, click Boundaries and choose Temperature.
2 Select Boundaries 3 and 5 only.
3 In the Settings window for Temperature, locate the Temperature section.
4 In the T0 text field, type Tc.

Temperature 3
1 In the Physics toolbar, click Boundaries and choose Temperature.
2 Select Boundary 1 only.
3 In the Settings window for Temperature, locate the Temperature section.
4 In the T0 text field, type Th-(Th-Tc)*s.

Temperature 4
1 In the Physics toolbar, click Boundaries and choose Temperature.
2 Select Boundary 4 only.
3 In the Settings window for Temperature, locate the Temperature section.
4 In the T0 text field, type Tc-(Tc-Th)*s.

MESH 1
Use a finer mesh setting to resolve the convection pattern well.

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 Finer.

9 | FREE CONVECTION IN A POROUS MEDIUM

STUDY 1

Step 1: Stationary
Set up an auxiliary continuation sweep for the alphap parameter to improve stability.

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 Click Add.
5 In the table, enter the following settings:

Parameter name Parameter value list Parameter unit

alphap (Fluid volumetric thermal 0 1e-12 1e-11 1e-10 1e-9 1/K
expansion) 1e-8 1e-7 6e-7 8e-7 1e-6

6 In the Home toolbar, click Compute.

RESULTS

Velocity (br)
The first default plot group shows the velocity magnitude. Refine the resolution for the
surface plot to get a smooth velocity field. Add an arrow plot to see the flow direction and
compare with Figure 3.

Surface
1 In the Model Builder window, expand the Velocity (br) node, then click Surface.
2 In the Settings window for Surface, click to expand the Quality section.
3 From the Resolution list, choose Finer.
4 In the Velocity (br) toolbar, click Plot.

Arrow Surface 1
1 In the Model Builder window, right-click Velocity (br) and choose Arrow Surface.
2 In the Velocity (br) toolbar, click Plot.

Pressure (br)
Because the cavity is closed, the pressure distribution is solely due to gravity.

Temperature (ht)
The third default plot group shows the temperature field as a surface plot (Figure 2).

10 | FREE CONVECTION IN A POROUS MEDIUM

Isothermal Contours (ht)
The fourth default plot group shows the temperature contours.

11 | FREE CONVECTION IN A POROUS MEDIUM

12 | FREE CONVECTION IN A POROUS MEDIUM