INTRODUCTION TO

Porous Media Flow Module

Introduction to the Porous Media Flow Module
© 1998–2019 COMSOL
Protected by patents listed on [Link]/patents, and U.S. Patents 7,519,518; 7,596,474; 7,623,991; 8,457,932;
8,954,302; 9,098,106; 9,146,652; 9,323,503; 9,372,673; 9,454,625; and 10,019,544. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement ([Link]/comsol-license-agreement) and may be used or copied only under the terms of the
license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Compiler, COMSOL Server,
and LiveLink are either registered trademarks or trademarks of COMSOL AB. All other trademarks are the property
of their respective owners, and COMSOL AB and its subsidiaries and products are not affiliated with, endorsed by,
sponsored by, or supported by those trademark owners. For a list of such trademark owners, see [Link]/
trademarks.
Version: COMSOL 5.5

Contact Information
Visit the Contact COMSOL page at [Link]/contact to submit general inquiries, contact
Technical Support, or search for an address and phone number. You can also visit the Worldwide
Sales Offices page at [Link]/contact/offices for address and contact information.

If you need to contact Support, an online request form is located at the COMSOL Access page at
[Link]/support/case. Other useful links include:

• Support Center: [Link]/support

• Product Download: [Link]/product-download
• Product Updates: [Link]/support/updates
• COMSOL Blog: [Link]/blogs
• Discussion Forum: [Link]/community
• Events: [Link]/events
• COMSOL Video Gallery: [Link]/video
• Support Knowledge Base: [Link]/support/knowledgebase

Part number: CM024802

Contents

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Porous Media Flow Module Interfaces . . . . . . . . . . . . . . . . . . . . . 6
Chemical Reaction and Mass Transport. . . . . . . . . . . . . . . . . . . . . . . 8
Fluid Flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
Porous Media Flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
Heat Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Structural Mechanics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Physics Interface Guide by Space Dimension and Study Type . . . 12
Tutorial Example — Packed Bed Latent Heat Storage. . . . . . . 16

|3
4|
Introduction

Porous media are encountered in many natural and man-made systems. The need
for advanced porous media modeling spans many industries and application areas
such as processes in fuel cells, drying of pulp and paper, food production, filtration
processes, and so on.
The Porous Media Flow Module extends the COMSOL Multiphysics modeling
environment to the quantitative investigation of mass, momentum and energy
transport in porous media. It is designed for researchers, engineers, teachers, and
students, and it suits both single-physics and multiphysics modeling.
The contents of the Porous Media Flow Module are a set of fundamental building
blocks which cover a wide array of physics questions. The physics interfaces it
offers work on their own or linked to each other. They can also be coupled to
physics interfaces already built into COMSOL Multiphysics, or to new equations
you create.
The physics interfaces, options, and functionalities in this module are tailored to
account for processes in porous media. The Heat Transfer interfaces, for example,
include options to automate the calculation of effective thermal properties for
multicomponent systems. The fluid flow equations represent a wide range of
possibilities. Included are Richards’ equation, which describes nonlinear flow in
variably saturated porous media. The options for saturated porous media include
Darcy’s law for slow flow and the Brinkman equations where shear is non
negligible. The Laminar Flow and Creeping Flow interfaces cover free flows at
different Reynolds numbers. The module also treats the transport of chemical
species. The Transport of Diluted Species interface account for the transport of
species in free flow, saturated, and partially saturated porous media. A number of
examples link these physics interfaces together.
You can use the Porous Media Flow Module in combination with essentially any
of the other products in the COMSOL Multiphysics product suite. The module is
a general interface for computing optimal solutions to engineering problems to,
for example, improve a design so that it minimizes energy consumption or
maximizes the output.

|5
Porous Media Flow Module Interfaces

The Porous Media Flow Module contains a number of physics interfaces that
predefine equations or sets of equations adapted to mass, momentum and energy
transport in porous media. You can take the equations in these physics interfaces
and their variables and then modify them, link them together, and couple them to
physics interfaces represented elsewhere in COMSOL Multiphysics.
The physics interfaces are based on the laws for conservation of mass, momentum,
and energy in porous media. The flow models contain different combinations and
formulations of the conservation laws that apply to the physics of the fields. These
laws of physics are translated into partial differential equations and are solved
together with the specified initial and boundary conditions.
A physics interface defines a number of features. These features are used to specify
the fluid properties, initial conditions, boundary conditions, and possible
constraints. Each feature represents an operation describing a term or condition in
the conservation equations. Such a term or condition can be defined on a
geometric entity of the component, such as a domain, boundary, edge (for 3D
components), or point.
Another useful tool in these physics interfaces is the ability to describe material
properties such as density and viscosity by entering expressions that describe them
as a function of other parameters, such as species concentration, pressure, or
temperature. Many materials in the material libraries use temperature- and
pressure-dependent property values.
Figure 1 shows the group of physics interfaces available with this module in
addition to the COMSOL Multiphysics basic license. Use these interfaces to
model chemical species transport, fluid flow, heat transfer and solid mechanics, to

6|
make modeling easier, something that is briefly discussed next. See also Physics
Interface Guide by Space Dimension and Study Type.

Figure 1: The physics interfaces for the Porous Media Flow Module as shown in the Model Wizard for a
3D application.

The module handles time-dependent and stationary problems for 1D, 2D, and 3D
geometries with axisymmetry for 1D and 2D. The predefined physics interfaces
cover four main categories: Chemical Species Transport ( ), Fluid Flow ( ),
Heat Transfer ( ), and Structural Mechanics ( ), as discussed in the next
pages.

|7
Chemical Reaction and Mass Transport
The Transport of Diluted Species interface ( ) simulates chemical species
transport through convection (when coupled to fluid flow), diffusion, and
reactions, for mixtures where one component, a solvent, is present in excess.
The Transport of Diluted Species in Porous Media interface ( ) is tailored to
model solute transport in saturated and partially saturated porous media. This
physics interface characterizes the rate and transport of individual or multiple and
interacting chemical species for systems containing fluids, solids, and gases. The
equations supply predefined options to describe mass transfer by convection,
adsorption, dispersion, diffusion, volatilization and reactions. You define the
convective velocity from either of the included physics interfaces, or you set it to
a predefined velocity profile.
The Laminar Flow, Diluted Species interface ( ) under the Reacting Flow
branch combines the functionality of the Single-Phase Flow and Transport of
Diluted Species interfaces. This multiphysics interface is primarily applied to
model flow at low to intermediate Reynolds numbers in situations where the mass
transport and flow fields are coupled.
The Transport of Diluted Species interface ( ) under the Reacting Flow in
Porous Media branch combines the Brinkman Equations and the Transport of
Diluted Species in Porous Media interfaces. This multiphysics interface is primarily
applied to model the transport of diluted reacting mixtures in porous media.
The Transport of Diluted Species in Fractures interface ( ) is used to model the
transport of a solutes along thin porous fractures, taking into account diffusion,
dispersion, convection, and chemical reactions. The fractures are defined by
boundaries in 2D and 3D, and the solutes are diluted in a solvent. The mass
transport equation solved along the fractures is the tangential differential form of
the convection-diffusion-reaction equation. Different effective diffusivity models
are available.
The Moisture Transport in Building Materials interface ( ) is used to compute
the relative humidity field in building materials. It simulates moisture transport by
taking in account moisture storage, liquid transport by capillary suction forces.
The Moisture Transport in Air interface ( ) is used to compute the relative
humidity distribution in air. It simulates moisture transport by vapor convection
and diffusion in moist air and the evaporation or condensation on walls.
The Laminar Flow interface under the Moisture Flow folder ( ) combines all
features from the Moisture Transport in Air and Single-Phase Flow interfaces. The
Moisture Flow multiphysics coupling is automatically added. The fluid properties
may depend on vapor concentration. The physics interface supports low Mach
numbers (typically less than 0.3). Similarly to Conjugate Heat Transfer coupling,

8|
the moisture transport interface can be coupled with additional turbulence models
for the free flow with the use of the CFD Module.

Fluid Flow
Flow in porous media normally occur at low Reynolds numbers. The Reynolds
number (Re) is a measure of the ratio of the fluid viscous to the inertial forces
acting on the fluid and is given by: Re=ρUL/μ, where ρ is the fluid density, U is
a characteristic velocity, L is a characteristic length scale, and μ is the dynamic
viscosity.
The Creeping Flow interface ( ) approximates the Navier-Stokes equations for
the case when the Reynolds number is significantly less than 1. This is often
referred to as Stokes flow and it is appropriate for use when viscous flow is
dominant.
The Laminar Flow interface ( ) describes fluid motion when turbulences are not
present, and the Reynolds number is less than approximately 1000. The physics
interface solves the Navier-Stokes equations for incompressible, weakly
compressible, or compressible flows where the Mach number (Ma) is less than 0.3.
The Phase Transport interface ( )under the Multiphase Flow branch is used to
simulate the transport of multiple immiscible phases in free flow. This interface
solves for the averaged volume fractions of the phases.

Porous Media Flow

The Darcy’s Law interface ( ) describes fluid movement through interstices in
a porous medium. This physics interface can be used to model low velocity flows,
for which the pressure gradient is the major driving force, and the flow is mostly
influenced by the frictional resistance within the pores. Its use is within very low
flows, or media where the permeability and porosity are very small. You can also
set up multiple Darcy’s Law interfaces to model multiphase flows involving more
than one mobile phase. The Darcy’s Law interface can also add non-Darcian
effects such as Ergun or Forchheimer drag terms, which are important for
Reynolds numbers bigger than 100.
The Fracture Flow interface ( ) is a variant of Darcy’s law that defines the flow
along interior boundaries representing fractures within a porous (or solid)
medium.
The Richards’ Equation interface ( ) analyzes flow in variably saturated porous
media. With variably saturated flow, the permeability changes as fluids move
through the medium, filling some pores and draining others. Richards’ equation
appears similar to the saturated flow equation set out in Darcy’s Law, but it is
notoriously nonlinear due to changes in material properties from unsaturated to

|9
saturated conditions. The analytic formulas of van Genuchten and Brooks and
Corey are frequently employed with variably saturated flow modeling.
The Brinkman Equations interface ( ) is used to model incompressible, weakly
compressible, or compressible flows where the Mach number (Ma) is less than 0.3.
You can select the Stokes-Brinkman flow feature to reduce the equations’
dependence on inertial effects, when the Reynolds number is significantly less than
one. The Brinkman Equations interface extends Darcy’s law to describe the
dissipation of the kinetic energy by viscous shear, similar to the Navier-Stokes
equation. The Brinkman Equations interface can also add a Ergun or Forchheimer
drag term, which are viscous drags proportional to the square of the velocity.
These non-Darcian effects are important for Reynolds numbers bigger than 100.
The Free and Porous Media Flow interface ( ) is useful for modeling problems
where free flow is connected to porous media, such as in fixed-bed reactors and
catalytic converters. The Free and Porous Media Flow interface is used over at least
two different domains, a free channel and a porous medium. The physics interface
adds functionality that allows the equations to be optimized according to the flow
properties of the relevant domain. For example, you can select the
Stokes-Brinkman flow feature to reduce the equations’ dependence on inertial
effects in the porous domain, or just the Stokes’ flow feature to reduce the
equations’ dependence on inertial effects in the free channel.
The Phase Transport in Porous Media interface ( ) is used to simulate the
transport of multiple immiscible phases through a porous medium. The interface
solves for the averaged volume fractions (saturations) of the phases, and it does not
track the interface between the different phases, although microscopic interfacial
effects are taken into account in the macroscopic equations through capillary
pressure functions.
The Multiphase Flow in Porous Media interface ( ) combines the functionality
of the Darcy’s Law and Phase Transport in Porous Media interfaces. This
multiphysics interface is intended to model flow and transport of multiple
immiscible phases in a porous medium.
As always, the physics interfaces give you direct access to defining, with either
constants or expressions, the material properties that describe the porous media
flow. This includes the density, dynamic viscosity, permeability, and porosity.

Heat Transfer
The Heat Transfer interfaces apply to systems consisting of solids, fluids, and
fluid-solid mixtures, and are able to calculate effective properties for porous media
consisting of several fluids, gases and solid components, such as a rock formation
with different mineral proportions.

10 |
The Heat Transfer in Solids interface ( ) describes, by default, heat transfer by
conduction. The physics interface is also able to account for the heat flux due to
translation in solids as well as for solid deformation, including volume or surface
changes. In the case of irreversible thermally induced transformation, it accounts
for the enthalpy and material properties changes.
The Heat Transfer in Fluids interface ( ) accounts for conduction and
convection in gases and liquids as the default heat transfer mechanisms. The
coupling to the flow field in the convection term may be entered manually, or it
may be selected from a list that couples the heat transfer to an existing Fluid Flow
interface. The Heat Transfer in Fluids interface can be solved simultaneously with
the Laminar Flow interface, or when the flow field has already been calculated and
the heat transfer problem is solved afterward, typically for simulations of forced
convection.
The Heat Transfer in Solids and Fluids interface ( ) contains solids and fluids
domains by default. It is aimed to simplify the setup of models where capabilities
of Heat Transfer in Solids interface ( ) and Heat Transfer in Fluids interface
( ) are used, in particular in conjugate heat transfer applications.
The Conjugate Heat Transfer branch ( ) combine features from the Heat
Transfer and Single-Phase Flow interfaces to describe heat transfer in solids and
fluids and non isothermal flow in fluids. The heat transfer process is tightly
coupled with the fluid flow problem via a predefined multiphysics coupling.
Additional turbulence models for the free flow are available with the use of the
CFD Module.
The Heat Transfer in Fractures interface ( ) under the Thin Structures branch
( ) describes heat transfer in fractures and thin porous media. It provides
efficient models defined at the boundaries level representing thin
three-dimensional domains. The simplest model assumes that the temperature
changes through the fracture thickness can be neglected, while the general model
computes the temperature variation across the fracture.
The Building Materials interface ( ) combines the Heat Transfer in Building
interface ( ) with the Moisture Transport in Building Materials interface ( ).
It can be used to model different moisture variations phenomena in building
components such as drying of initial construction moisture, condensation due to
migration of moisture from outside to inside, or moisture accumulation by
interstitial condensation due to diffusion.
The Moist Air interface ( ) combines the Heat Transfer in Moist Air interface
( ) with the Moisture Transport in Air interface ( ). It is used to simulate the
coupling between heat transfer and vapor transport in air and the evaporation and
condensation on walls.

| 11
The Heat Transfer in Porous Media interface ( ) combines the heat conduction
and convection in a solid-fluid system. This physics interface provides mixing rules
for calculating the effective heat transfer properties, expressions for heat dispersion
in porous media. Dispersion is caused by the tortuous path of the liquid in the
porous medium, which would not be described if only the mean convective term
was taken into account. This physics interface may be used for a wide range of
porous materials, from porous structures to the simulation of heat transfer in soils
and rocks, and also to model heat transfer in fractures.
The Local Thermal Nonequilibrium (LTNE) multiphysics interface ( ) is
intended to simulate heat transfer in porous media where the temperatures in the
porous matrix and the fluid are not in equilibrium. The absence of thermal
equilibrium can result from fast transient changes, but it can also be observed in
stationary scenarios. Typical applications are rapid heating (or cooling) of porous
media using a hot (or cold) fluid, or internal heat generation in one of the phases.
This effect can be observed in nuclear devices, electronics systems, and fuel cells.
These features interact seamlessly and can be used in combination in a single
application. Surface-to-surface radiation can also be included in the energy
equation, although this requires a license for the Heat Transfer Module.

Structural Mechanics
The Poroelasticity multiphysics interface ( ) combines a transient formulation
of Darcy’s law with a linear elastic material included in the Solid Mechanics
interface. The poroelasticity coupling means that the pore fluid affects the
compressibility of the porous medium, as well as changes in volumetric strains
affect the porosity and the fluid flow.

Physics Interface Guide by Space Dimension and Study Type

The table lists the physics interfaces available with this module in addition to those
included with the COMSOL basic license.

PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE

DIMENSION

Chemical Species Transport

Transport of Diluted tds all dimensions stationary; time dependent

Species
Transport of Diluted tds all dimensions stationary; time dependent
Species in Porous Media

12 |
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION

Transport of Diluted dsf 3D, 2D, 2D stationary; time dependent

Species in Fractures axisymmetric
Moisture Transport

Moisture Transport in mt all dimensions stationary; time dependent

Building Materials
Moisture Transport in mt all dimensions stationary; time dependent
Air
Moisture Flow

Laminar Flow — 3D, 2D, 2D stationary; time dependent

axisymmetric
Reacting Flow

Laminar Flow, Diluted — 3D, 2D, 2D stationary; time dependent

Species axisymmetric
Reacting Flow in Porous Media

Transport of Diluted rfds 3D, 2D, 2D stationary; time dependent

Species axisymmetric
Fluid Flow

Single-Phase Flow

Creeping Flow spf 3D, 2D, 2D stationary; time dependent

axisymmetric
Laminar Flow spf 3D, 2D, 2D stationary; time dependent
axisymmetric
Phase Transport

Phase Transport phtr 3D, 2D, 2D stationary; time dependent

axisymmetric

| 13
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION

Porous Media and Subsurface Flow

Brinkman Equations br 3D, 2D, 2D stationary; time dependent

axisymmetric
Darcy’s Law dl all dimensions stationary; time dependent

Fracture Flow esff 3D, 2D, 2D stationary; time dependent

axisymmetric
Richards’ Equation dl all dimensions stationary; time dependent

Multiphase Flow in — 3D, 2D, 2D stationary; time dependent

Porous Media axisymmetric
Free and Porous Media fp 3D, 2D, 2D stationary; time dependent
Flow axisymmetric
Phase Transport in phtr 3D, 2D, 2D stationary; time dependent
Porous Media axisymmetric
Heat Transfer

Heat Transfer in Solids ht all dimensions stationary; time dependent

Heat Transfer in Fluids ht all dimensions stationary; time dependent

Heat Transfer in Solids and ht all dimensions stationary; time dependent

Fluids
Thin Structures

Heat Transfer in htlsh 3D, 2D, 2D stationary; time dependent;

Fractures axisymmetric thermal perturbation,
frequency domain
Heat and Moisture Transport

Building Materials — all dimensions stationary; time dependent;

thermal perturbation,
frequency domain

14 |
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION

Moist Air — all dimensions stationary; time dependent;

thermal perturbation,
frequency domain
Heat and Moisture Flow

Laminar Flow — 3D, 2D, 2D stationary; time dependent

axisymmetric
Heat Transfer in Porous ht all dimensions stationary; time dependent
Media
Local Thermal — all dimensions stationary; time dependent;
Nonequilibrium thermal perturbation,
frequency domain
Structural Mechanics

Poroelasticity poro 3D, 2D, 2D stationary; time dependent

axisymmetric
1 This physics interface is included with the core COMSOL package but has added functionality
for this module.
1
Requires the addition of the Heat Transfer Module or CFD Module.
2 Requires the addition of the CFD Module.

| 15
Tutorial Example — Packed Bed Latent Heat Storage

Thermal energy storage (TES) units are used to accumulate thermal energy from
solar, geothermal, or waste heat sources. The simplest TES units are built from
water tanks, where the solar energy is stored as sensible heat. These systems are
called sensible heat storage (SHS) units. The thermal capacity of these tanks can
be further increased by including latent heat, which gives rise to latent heat storage
(LHS) units. Typically, LHS tanks contain spherical capsules filled with paraffin
wax as phase change material. Paraffin is a suitable phase change material, as it is
relatively inexpensive, reliable and nontoxic, and it is commercially available for a
wide rage melting temperatures.

This example is inspired by the experimental investigation found in Ref. 1. It

models the flow through a packed-bed storage tank, and it includes the effects of
heat transfer with phase change and local thermal nonequilibrium (LTNE) while
charging the LHS unit.
Inlet

Glass wool
insulation

Spherical capsules
filled with paraffin
Outlet

Figure 2: Latent heat storage unit.

16 |
M ODEL DEFINITION
Paraffin-filled spherical capsules with a diameter of dp = 55 mm are stored in a tank
of 36 cm in diameter and 47 cm in height. The porosity of this bed is εp = 0.49.
The model geometry is shown in Figure 2.
The material properties of paraffin wax are listed in the following table.

MATERIAL PROPERTY PARAFFIN, SOLID PARAFFIN, LIQUID

Density ρ (kg/m3) 861 778

Heat Capacity, Cp (J/(kg·K)) 1850 2384
Thermal conductivity, k (W/(m·K)) 0.4 0.15
Melting temperature, Tm (°C) 60
Latent heat of fusion, L (J/kg) 213

Geometry, material properties, and operating conditions are taken from Ref. 1.
The initial temperature in the tank is set to 32 °C. Warm water flows through the
tank with a flow rate of min = 2 l/min. During thermal charging the water is
continuously heated up by a solar collector that delivers a power of Qu = 375 W.
The temperature difference at the tank’s inlet and outlet is given by the relation
Qu
--------
- = ρC p ( T in – T out ) (1)
m in

here, Tin and Tout are the inlet and outlet temperatures, and ρ and Cp are the
density and heat capacity of water.
Ergun equation describes non-Darcian flow through the packed bed, which
estimates the pressure drop as a function of the velocity field u
μ 1,75 ( 1 – ε p )
∇p = – --- u – ------------------------------
-ρ u u (2)
κ dp εp
3

here, μ (Pa·s) and ρ (kg/m3) are the viscosity and density of water, dp (m) is the
spheres’ diameter, and εp the bed porosity. The permeability κ (m2) of the packed
bed is given by
2 3
dp εp
κ = --------------------------------2 (3)
150 ( 1 – ε p )

The Reynolds number in the packed bed can be estimated as

d p vρ
Re = ----------------------- (4)
( 1 – ε p )μ

| 17
The maximum velocity in the tank, v, is about 6 mm/s, which implies a Reynolds
number of about 600. For this Reynolds number the flow field is assumed to be
independent of the temperature distribution, such that a stationary field can be
computed before running the thermal simulation. This is a reasonable
simplification that reduces the computational effort.
The relative large diameter of the capsules as compared to the tank dimensions
suggests a significant temperature difference between the encapsulated paraffin
and the surrounding water flow, thus a local thermal nonequilibrium (LTNE)
approach is considered in the example.
The Local Thermal Nonequilibrium multiphysics interface is used to couple the
Heat Transfer in Solids with the Heat Transfer in Fluids interfaces. The heat
transferred from the paraffin-filled capsules to the water is modeled with a heat
source
q sf
Q f = ------ ( T s – T f ) (5)
εp

here, Ts and Tf are the paraffin and water temperatures, and qsf (W/(m3·K)) is the
interstitial convective heat transfer coefficient, which for spherical capsules reads
6 ( 1 – εp )
q sf = ----------------------- h sf (6)
dp

The interstitial heat transfer coefficient hsf follows a Nusselt number correlation
see the Theory for the Local Thermal Nonequilibrium Interface section in the
documentation for more information. Convection inside the capsules is neglected,
thus the paraffin wax is treated as a solid or immobile liquid.

R ESULTS
The tank reaches a temperature of 70 °C after approximately 13 hours. The
resulting velocity and temperature distribution is shown in Figure 3

18 |
Figure 3: Velocity field (streamlines) with the gray color indicating the pressure and temperature field
(color) after 12 hours of thermal charging.

Figure 4 shows the evolution of the paraffin temperature, water temperature and
the weighted average (porous medium) temperature at three different points
located in the central axis of the tank. During the phase change, the encapsulated
paraffin is not in thermal equilibrium with the surrounding water. Measuring the
water temperature alone at the inlet or outlet does not give accurate information
about neither the temperature inside the capsules nor the phase in which the
paraffin wax is.

| 19
Figure 4: Evolution of water (dashed), paraffin (dotted) and average porous medium temperature (solid)
during phase change for top (red), center (green) and low (blue) position.

The evolution of the phase distribution changes as the paraffin wax warms up.
Liquid paraffin appears after approximately 4 hours, when water heats up the
paraffin wax to its melting temperature of 60°C. The wax is completely molten
after about 10 hours. The latent heat storage tank is considered fully charged as
soon as a temperature of 70°C is reached everywhere in the tank, which happens
after approximately 13 hours.

20 |
Figure 5: Paraffin solid phase (blue) and liquid phase (yellow) after 7 hours of thermal charging the LHS
unit.

Figure 5 shows the paraffin phase distribution after 7 hours. Near the walls, where
the flow velocity is negligible, the phase transition has not yet begun and the
paraffin wax is still solid, while it is melted in the center of the tank.

R EFERENCES
1. N. Nallusamy et al, “Study on performance of a packed bed latent heat thermal energy
storage unit integrated with solar water heating system,” Journal of Zhejiang
University-SCIENCE A, vol. 7, pp. 1422–1430, 2006.

Model Wizard

Note: These instructions are for the user interface on Windows but apply, with
minor differences, also to Linux and Mac.
1 To start the software, double-click the COMSOL icon on the desktop. When
the software opens, you can choose to use the Model Wizard to create a new

| 21
 COMSOL model or Blank Model . to create one manually. For this tutorial,
click the Model Wizard button.
If COMSOL is already open, you can start the Model Wizard by selecting
New from the File menu and then click Model Wizard .
The Model Wizard guides you through the first steps of setting up a model. The
next window lets you select the dimension of the modeling space.
2 In the Select Space Dimension window click 2D Axisymmetric .
3 In the Select Physics tree, expand the Fluid Flow folder, and under the
Porous Media and Subsurface Flow folder double-click
Free and Porous Media Flow to add it to the Added physics interfaces list.
You can also click Add or right-click the icon and choose Add Physics .
4 In the Select Physics tree under the Heat Transfer folder, double-click
Local Thermal Nonequilibrium to add it to the Added physics interfaces
list. You can also click Add or right-click and choose Add Physics .
5 Click Done .

Geometry

You can build the geometry from geometric primitives. Here, use instead a file
containing the sequence of geometry features that has been provided for
convenience.
The location of the file used in this exercise varies based on your installation. For
example, if the installation is on your hard drive, the file path might be similar to
C:\Program Files\COMSOL\COMSOL55\Multiphysics\applications\.

1 On the Geometry toolbar choose Import .

2 In the Settings window for Import, locate the Import section, and browse to
the model’s Application Libraries folder. Double-click the file
packed_bed_latent_heat_storage.mphbin.
3 In the Settings window for Import, click Import.

22 |
4 Click Build All Objects to display the geometry.

Global Definitions — Parameters

Add the parameters needed in this example.

Parameters
1 In the Model Builder window under Global Definitions click
Parameters .
2 In the Settings window for Parameters enter the following data

| 23
Mate rials

Define the material nodes for water, paraffin, and the insulating glass.

Water, liquid
1 In the Home toolbar, click
Add Material .
2 Go to the Add Material window,
expand the Built-in folder and select
Water, liquid. Click
Add to Component in the window
toolbar to add it o the Materials list.
You can also double-click on it or
right-click and choose
Add to Component 1 (comp 1) .
3 On the Home toolbar click
Add Material again to close the
window.

Paraffin, solid
4 In the Model Builder window, under
Component 1 (comp 1) ,
right-click Materials and choose
Blank Material .
5 In the Label text field under the Settings window for Material, type
Paraffin, solid.

Repeat the same procedure to create two more material nodes, call them
Paraffin, liquid and Glass_wool respectively. After setting up the physics
interfaces you can fill in the required material properties in these nodes.

Definitions

Average 1 (aveop1)
Create an average operator to evaluate the inlet temperature.
1 In the Definitions toolbar click Nonlocal Couplings and choose
Average .

24 |
2 In the Settings window for Average locate the Source Selection section, and
from the Geometric entity level list choose Boundary.
3 Select Boundary 2 only.

Minimum 1 (minop1)
Create a minimum operator to evaluate the minimum temperature in the packed
bed.
4 In the Definitions toolbar click Nonlocal Couplings and choose
Minimum .
5 Select Domain 2 only.

Variables 1
As water is pumped through a closed loop, its temperature increases over time
during the charging process. Define the variables deltaT and T_in as functions
of the outlet temperature using Equation 1. Also, define the minimum temperature
in the packed bed with the help of the minimum operator minop1, which is used
later as stop condition in the solver settings.
6 In the Definitions toolbar click Local Variables .
7 In the Settings window for Variables, locate the Variables section and enter the
following settings.

The variables [Link] and [Link] are the heat capacity and density of water as
defined by the Heat Transfer in Fluids interface.

Explicit 1
Create a selection for the tank outer boundary to apply a heat flux condition.

| 25
8 In the Definitions toolbar click Explicit .
9 In the Settings window for Explicit locate the Input Entities section, and from
the Geometric entity level list choose Boundary.
10Select the Group by continuous tangent check box.
11Select Boundaries 14–17, 21–23, 26, 27, 30, and 31. In the Label text field type
Heat flux boundary.

Free and Porous Media Flow

Fluid and Matrix Properties1

1 In the Model Builder window, under Component 1 (comp 1) , click
Free and Porous Media Flow (fp) and select Domains 1–3. These are the
free flow and packed bed domains.
2 In the Physics toolbar click Domains and choose
Fluid and Matrix Properties . Select Domain 2 only.
3 In the Settings window for Fluid and Matrix Properties locate the
Fluid Properties section. From the Fluid material list choose
Water, liquid (mat1).
4 Locate the Porous Matrix Properties section. From the εp list choose
User defined and type por.
5 From the Permeability model list choose Non-Darcian, and from the
Non-Darcian flow model list choose Ergun. In the Particle diameter dp field
type dp.
Next, apply the boundary conditions for the flow.

Inlet 1
1 In the Physics toolbar click Boundaries and choose Inlet . Select
Boundary 7 only.
2 In the Settings window for Inlet locate the Boundary Condition section, from
the list choose Fully developed flow. Click the Flow rate button, and in the
Flow rate V0 field type m_in.

Outlet 1
1 In the Physics toolbar click Boundaries and choose Outlet . Select
Boundary 2 only. The default outlet condition specifies a zero relative pressure,
and suppresses back flow effects.

26 |
Heat Transfer in Solids

Solid 1
Continue with setting up the heat transfer properties for the insulating glass wool.
1 In the Model Builder window, under Component 1 (comp1) click
Heat Transfer in Solids (ht) and select Domains 2 and 4 only.
2 In the Model Builder window, under
Component 1 (comp1)>Heat Transfer in Solids (ht) click Solid 1 .
3 In the Settings window for Solid, locate the Solid Material section and from the
list, choose Glass wool (mat4).

Solid 2
Continue with setting up the paraffin properties.
4 In the Physics toolbar click Domains and choose Solid . Select Domain
2 only.
5 In the Physics toolbar click Attributes and choose
Phase Change Material .
6 In the Settings window for Phase Change Material locate the Density section.
7 From the ρ list, choose User defined. In the associated text field type 820[kg/
m^3], which is the average density of liquid and solid paraffin. The paraffin
density is assumed to remain constant, which is a reasonable assumption.
8 Locate the Phase Change section. In the Tpc, 1 → 2 field type 60[degC]. In the
in the ΔT1 → 2 field type 2[K], and in the L1 → 2 field type 213[kJ/kg].
9 Locate the Phase 1 section. From the Material, phase 1 list choose
Paraffin, solid (mat2).
10Locate the Phase 2 section. From the Material, phase 2 list choose
Paraffin, liquid (mat3).

Initial Values 1
11In the Model Builder window, click Initial Values 1. In the Settings window
locate the Initial Values section and type T0 in the T text field.

Heat Flux 1
Apply the boundary conditions for the heat transfer in the packed bed.
12In the Physics toolbar click Boundaries and choose Heat Flux .
13In the Settings window for Heat Flux locate the Boundary Selection section,
from the Selection list choose Heat flux boundary.

| 27
14Locate the Heat Flux section and click the Convective heat flux button. In the
Heat transfer coefficient h field type 5[W/m^2/K].

Heat Transfer in Fluids 2 (ht2)

Continue with setting up the heat transfer properties for the water flow.
1 In the Model Builder window under Component 1 (comp1) click
Heat Transfer in Fluids 2 (ht2) and select Domains 1–3 only.

Fluid 1
2 In the Model Builder window, under
Component 1 (comp1)>Heat Transfer in Fluids 2 (ht2) click Fluid 1 .
3 In the Settings window for Fluid, locate the Fluid Material section and from the
list choose Water, liquid (mat1).

Initial Values 1
4 In the Model Builder window, under
Component 1 (comp1)>Heat Transfer in Fluids 2 (ht2) click
Initial Values 1 .
5 In the Settings window for Initial Values, locate the Initial Values section.
6 In the T2 text field type T0.
Next, apply the boundary conditions for heat transfer in water.

Outflow 1
7 In the Physics toolbar click Boundaries and choose Outflow . Select
Boundary 2 only.

Temperature 1
8 In the Physics toolbar click Boundaries and choose Temperature . Select
Boundary 7 only, which corresponds to the water inlet.
9 In the Settings window for Temperature locate the Temperature section, in the
T0 field type T_in.
Now, fill out the remaining material properties. Since you set up the physics, the
software automatically detects which properties are required for the simulation.

Paraffin, solid (mat2)

10In the Model Builder window, under Component 1 (comp1)>Materials
click Paraffin, solid (mat2).

28 |
11In the Settings window for Material locate the Material Contents section and
enter the following settings:

Paraffin, liquid (mat3)

12In the Model Builder window, under Component 1 (comp1)>Materials
click Paraffin, liquid (mat3).
13In the Settings window for Material locate the Material Contents section and
enter the following settings:

Glass wool (mat4)

14In the Model Builder window, under Component 1 (comp1)>Materials
click Glass wool (mat4) and select Domain 4 only.
15In the Settings window for Material locate the Material Contents section and
enter the following settings:

Add Multiphysics

1 In the Physics toolbar, click Add Multiphysics to couple the

Free and Porous Media Flow with the Heat Transfer in Fluids
interface.

| 29
2 Go to the Add Multiphysics window and find the
Select the physics interfaces you want to couple subsection. In the table, enter
the following settings:

3 In the tree, select Fluid Flow>Nonisothermal Flow>Laminar Flow .

4 Click Add to Component. In the Physics toolbar, click again
Add Multiphysics to close the window.

Local Thermal Nonequilibrium 1 (ltne1)

5 In the Model Builder window under Component 1 (comp1)>Multiphysics
click Local Thermal Nonequilibrium 1 (ltne1) .
6 In the Settings window for Local Thermal Nonequilibrium locate the
Local Thermal Nonequilibrium Settings section, and in the θp text field type
1-por.

30 |
7 From the Interstitial convective heat transfer coefficient list choose
Spherical pellet bed, and in the rp text field type dp/2.

Me sh 1

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,
and from the Element size list choose Fine.
3 Click Build All , to see the automatically generated mesh.

| 31
The figure below shows the boundary layer mesh at the walls. Zoom in to the
mesh using the zoom function on the Graphics toolbar to confirm that it
matches the figure.

Add Study

Stationary
Add a stationary study step to solve for the velocity profile.
1 In the Home toolbar click Add Study . In the Add Study window find the
Studies subsection, expand the General Studies folder, and select
Stationary .

32 |
2 Find the Physics interfaces in study subsection, and the
Multiphysics couplings in study subsection. Enter the following settings in the
tables to solve fist the stationary flow field:

3 Click Add Study in the window toolbar.

4 To clear the default plots generated by the solver, in the Model Builder window
click Study 1 .
5 In the Settings window for Study locate the Study Settings section, and clear the
Generate default plots check box.

Next, add a Time Dependent step to solve the coupled heat transfer problem.

Time Dependent
6 In the Study toolbar click Study Steps and choose Time Dependent .
7 In the Settings window for Time Dependent locate the Study Settings section,
and from the Time unit list choose h.
8 In the Times text field type
range(0,0.25,3.75) range(4,5[min],8.75) range(9,0.25,24). This

| 33
 will carry on the simulation time for 24 hours. The time steps are chosen such
that the phase change is properly resolved.
9 From the Tolerance list choose User controlled, and in the Relative tolerance
text field type 1e-4.
10Locate the Physics and Variables Selection section and enter the following
settings in the tables.

11In the Study toolbar, click Show Default Solver .

12In the Solver Configurations node, expand the Solution 1 (sol1) node
and click the Time-Dependent Solver 1 node.
13In the Settings window for Time-Dependent Solver, click to expand the
Time Stepping section, and from the Steps taken by solver list choose Strict.
This option forces the time-dependent solver to strictly use the time steps
specified before.
The time it takes to charge the tank is not known a priori. In order to reduce
the computational time, a stop condition for the time dependent solver is used
to stop the simulation after a temperature of 70°C is reached in every point
inside the packed bed.
14Right-click the Time-Dependent Solver 1 node and choose
Stop Condition .
15In the Settings window for Stop Condition, locate the Stop Expressions section
and click Add.

34 |
16In the table enter the following settings.

17In the Home toolbar, click Compute . A warning message appears and states
that the stop condition was fulfilled after about 47000 seconds (about 13
hours).

Results

As we removed the default plot generation from the solver settings, follow the
steps to generate Figure 3 and Figure 4.

2D P LOT G ROUP 1
1 In the Home toolbar, click Add Plot Group and choose
2D Plot Group . In the Settings window for 2D Plot Group, type
Temperature and velocity field in the Label text field.
2 Click to expand the Title section, and from the Title type list choose None.

Surface 1
3 Right-click Temperature and velocity field and choose Surface .
4 In the Settings window for Surface locate the Expression section, and in the
Expression text field type T.
5 Locate the Coloring and Style section. From the Color table list choose
ThermalLight.
6 Right-click Surface 1 and choose Selection . Select Domain 4 only. This
plots the temperature field in the container.

| 35
Surface 2
7 In the Model Builder window, right-click the
Temperature and velocity field node and choose Surface . You can also
click Surface in the Temperature and velocity field toolbar.
8 In the Settings window for Surface, locate the Expression section, and in the
Expression text field type ltne1.T. This plots for the average temperature in the
packed bed.
9 Click to expand the Inherit Style section, and from the Plot list choose
Surface 1.

Streamline 1
10Right-click the Temperature and velocity field node and choose
Streamline . You can also click Streamline in the
Temperature and velocity field toolbar.
11In the Settings window for Streamline, locate the Selection section and select
Boundary 7 only. This is the inlet on the upper boundary.
12Locate the Coloring and Style section. Find the Point style subsection, and
from the Type list choose Arrow.
13Right-click the Streamline 1 node and choose Color Expression .
14In the Settings window for Color Expression locate the Expression section, and
in the Expression text field type p.
15Locate the Coloring and Style section, and from the Color table list choose
GrayScale. Select the Reverse color table check box.

1D P LOT G ROUP 2
Add a new data set to create Figure 4. Then evaluate and plot the different
temperatures over time with the help of table graphs.

Cut Point 2D 1
1 In the Results toolbar, click Cut Point 2D .
2 In the Settings window for Cut Point 2D, locate the Point Data section. In the
r text field type 0, and in the z text field type 0.05 0.235 0.42. These points
will evaluate temperatures in the packed bed at three different positions along
the symmetry axis.
3 In the Results toolbar, click Point Evaluation .
4 In the Settings window for Point Evaluation locate the Data section, and from
the Dataset list choose Cut Point 2D 1.

36 |
5 Locate the Expressions section. In the table enter the following settings:

6 Click Evaluate . A table is created in the window toolbar with the data for
the temperature variables evaluated at the cut-points for all the time steps.

Table Graph 1
7 Go to the Table window and click Table Graph in the window toolbar. This
generates a 1D Plot Group 2 node with a Table Graph 1 node.
8 In the Settings window for Table Graph locate the Data section, and from the
Plot columns list choose Manual.
9 In the Columns list choose
Paraffin temperature (K), Point: (0, 0.05),
Paraffin temperature (K), Point: (0, 0.235), and
Paraffin temperature (K), Point: (0, 0.42).
10Locate the Coloring and Style section. Find the Line style subsection, and from
the Line list choose Dotted.

Table Graph 2
11Right-click the Table Graph 1 node and choose Duplicate .
12In the Settings window for Table Graph locate the Data section, and in the
Columns list choose Water temperature (K), Point: (0, 0.05),
Water temperature (K), Point: (0, 0.235), and
Water temperature (K), Point: (0, 0.42).
13Locate the Coloring and Style section. Find the Line style subsection, and from
the Line list choose Dashed. From the Color list choose Cycle (reset).

| 37
Table Graph 3
14Right-click Table Graph 2 node and choose Duplicate .
15In the Settings window for Table Graph locate the Data section, and in the
Columns list choose
Porous medium temperature (K), Point: (0, 0.05),
Porous medium temperature (K), Point: (0, 0.235), and
Porous medium temperature (K), Point: (0, 0.42).
16Locate the Coloring and Style section. Find the Line style subsection and from
the Line list choose Solid.

1D Plot Group 2
17In the Model Builder under Results click 1D Plot Group 2 .
18In the Settings window for 1D Plot Group, type Temperature evolution in
the Label text field to rename the plot.
19Locate the Axis section and select the Manual axis limits check box. In the
x minimum text field, type 3.5, and in the x maximum text field type 9. This
zooms into the time range between 3.5 and 9 hours. In the y minimum text
field type 328, and in the y maximum text field type 344. to adjust the
temperature range.
20Click Plot in the Temperature evolution toolbar to visualize the data.
You clearly see that paraffin and water are not in thermal equilibrium, especially
during phase change of paraffin.

38 |