PDF
Free Convection in a Porous Medium
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.
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
(1)
.
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
(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 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.
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.
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
Pr
mu*Cp/kth
0.7
Prandtl number
Ra
Cp*rho^2*g_const*alphap*(Th-Tc)*L^3/(kth*mu)
1.5102E5
Rayleigh number
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.
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.
Fluid 1
1
In the Model Builder window, under Component 1 (comp1)>Heat Transfer in Porous Media (ht)>Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, 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 kf list, choose User defined. In the associated text field, type kth.
6
Locate the Thermodynamics, Fluid section. From the ρf list, choose User defined. In the associated text field, type rho.
7
From the Cp,f 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.
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type epsilon.
4
Locate the Heat Conduction, Porous Matrix section. From the kb list, choose User defined. Locate the Thermodynamics, Porous Matrix section. From the ρb list, choose User defined. From the Cp,b list, choose User defined.
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.
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 expansion)
0 1e-12 1e-11 1e-10 1e-9 1e-8 1e-7 6e-7 8e-7 1e-6
1/K
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).
Isothermal Contours (ht)
The fourth default plot group shows the temperature contours.