Excerpt from the Proceedings of the COMSOL Conference 2008 Hannover
Comsol Multiphysics in Education – Chemical Reactions, Heat
and Mass Transfer
R. Geike
TFH Berlin, Maschinenbau, Verfahrens- und Umwelttechnik
Corresponding author: TFH Berlin, FB VIII, Luxemburger Str. 10, D-13353 Berlin,
[Link]@[Link]
Abstract: In our Master program entitled - Process Simulation,
“Verfahrenstechnik / Process Engineering” the - Transport Processes with Comsol Multi-
main focus is placed on the use of computational physics.
methods, including computer exercises. This The education in these courses is
paper will provide the plan of a particular course complemented by the course entitled “Numerics
entitled “Transportprozesse” (Transport Proces- and Optimization”.
ses or Transport Phenomena). For the course entitled “Transport Processes”
The contents of the lecture and exercises its contents, organization and aim as well as the
encompass mass and heat transfer, chemical priorities of the exercises and the assessment of
reactions and fluid mechanics. That includes the students’ performance are all described in the
“single physic” and “multi physic”, steady state following sections - with the aid of three
and transient processes. The selected problems examples.
are solved analytically and numerically. The
software called Comsol Multiphysics is used to 2. The Course entitled “Transport Pro-
find the numerical solution. The organization of cesses”
the course, its aims, the priorities of the exercises
and the assessment of the students’ performance This course consists of 4 hours of instruction
are explained with the aid of three examples. per week: 2 hours of lecture and 2 hours of
computer exercises. Topics are taken from the
Keywords: mass and heat transfer, chemical fields of heat and mass transfer, chemical
reactions, fluid mechanics, Comsol Multi- reaction engineering and fluid dynamics:
physics, education. - steady state heat conduction - in a fin for
increasing heat emission and the
1. Introduction accompanying heating of a pipe,
- isothermal pore diffusion to obtain the effec-
Processes in Chemical Engineering and in tiveness factor of a catalyst particle with
other engineering disciplines are described by various shapes,
differential equations. Analytical solutions do - non-isothermal pore diffusion with a combi-
exist only for special cases, for simple nation of mass and heat transfer,
geometries and idealized conditions. The modern - transient processes of heating / cooling, e.g.
engineer needs numerical solutions for process the heating of a sphere and a one-sided fired
and apparatus design. wall or freezing of a piece of meat with a
Present day education of engineers has to phase change,
teach on the one hand the basics of various - axial dispersion in chemical reactors,
disciplines including traditional solution methods combined with the estimation of conversion,
and on the other hand the use of modern - electrically heated metallic component com-
numerical solution methods. bined with the heat emission,
Based on the course entitled “The Finite - laminar flow of Newtonian fluid to compare
Element Method” in the Bachelor program we go with Hagen-Poiseuille and non-Newtonian
more deeply into this field in our Master fluid in a tube,
program “Verfahrenstechnik / Process - flow and heat transfer over a step.
Engineering”. Three engineering courses deal In particular the basic and the chemical
with computational methods and include engineering modules of Comsol Multiphysics are
computer exercises: employed.
- Computational Fluid Dynamics with CFX,
Three points are used for the assessment of data tables for properties of materials, the change
the students’ performance. The first part consists of phase as well as other aspects.
of the grades for the reports of about 10 4) Alternatives are discussed for description
computer exercises. The students have to finish of problems: use of the axial dispersion model
the calculations and present and discuss the instead of the solution of the flow problem or the
results in a short paper. The second part of the use of different geometries (1D, 2D, 3D, use of
assessment is the grade for the written final symmetry).
exam, which focuses on the basics of transport 5) One aim of learning in this course
phenomena and of numerical solutions. The final concerns finding the similarities between various
part of the assessment is a grade for the processes. For example, heat transfer in a fin
individual homework assignment presented in without a source term is described by the same
the form of a paper and a Power Point differential equation that applies to a first order
presentation. Each student selects a special chemical reaction in a catalyst plate.
problem from a provided list, collects data and 6) Different results are obtained through
solves the problem on his own. variation of the number of elements. There is
discussion on the accuracy of calculations
3. Aim of Learning / Methodical Concept depending on the number of elements. The
dependence on element size can be compared
The methodical concept of the course with the influence of element size in the Euler or
concentrates on the following priorities: Runge-Kutta method.
1) At first very simple examples are solved
with basic geometry under idealized conditions. 4. Examples
This allows for analytical solutions - with pencil 4.1 Overview
and paper or with the help math and engineering
software entitled MAPLE. These analytical Two kinds of examples are employed: very
solutions can be compared with the results of simple examples for comparison with analytical
FEM calculations (Comsol Multiphysics) and solutions and slightly more complicated
with the results of simple numerical methods examples. The simple examples are concerned
(EULER method with EXCEL). with heat conduction in a fin, with isothermal
All possibilities are employed for the pore diffusion in a plate or a sphere, with axial
verification of the numerical results - e.g. the dispersion in chemical reactors, with laminar
comparison between the subdomain integral over flow of Newtonian fluid in a tube or with the
the reaction rate and the boundary integral over transient processes of heating / cooling of simple
the mass or heat flux. geometries.
2) Based on the simple examples students The second kind of examples uses slightly
progress to slightly more complicated geometries more complicated geometries or conditions.
or conditions. On the one hand, that shows the These examples might be the completion of the
path to the solution of practical problems. On the listed simple examples - accompanying heating
other hand it allows students to see the influence of a pipe, pore diffusion in different geometries
of simplifications / idealizations and to assess or the freezing of meat with phase change.
approximation formulas from literature. Such Examples are also used from the Comsol
simplifications may concern the isolation of parts Multiphysics model library, e. g. the flow with
of surface, the plug flow in reactors or heat heat transfer over a step or the electrically heated
exchangers or the flow behavior of a Bingham metallic component.
substance.
Interestingly the simplifications necessary for 4.2 Steady State Heat Conduction
analytical solutions may be a serious barrier for
the numerical solution with finite element The first example deals with heat conduction
methods. in a fin for increasing heat emission. The
3) Each example is used to teach new aspects simplified energy balance takes into account the
of modeling in Comsol Multiphysics in parti- heat conduction in the direction of longitudinal
cular regarding steady state and transient axis and heat emission from the surface of the
calculations, the import of geometry, the use of
fin. This results in the following ordinary The next step concerns optimization that
differential equation means the minimization of heat consumption for
T ′′ − m 2 ⋅ T = 0 maintaining the required product temperature.
and as a solution
cosh(m( L − x))
T ( x) = To
cosh(mL)
Through differentiation the following heat
flux at position x is obtained from this equation
for temperature:
sinh(m( L − x))
Q& ( x) = λ ⋅ A ⋅ To ⋅ m ⋅
cosh(mL)
Q& all = Q& ( x = 0) = λ ⋅ A ⋅ To ⋅ m ⋅ tanh(mL)
For x = 0 it is the heat flux which is fed to the
fin. Figure 1 shows the temperature and heat flux
Figure 2. Temperature field of the accompanying
for one example.
heating of a pipe.
105 300
4.3 Diffusion in Heterogeneous Catalysts
100 250
temperature [°C]
heat flux [W]
95 200 For the first order reaction in a heterogeneous
90 150 catalyst plate with diffusion transport and
85 100
consumption by a chemical reaction the same
T mathematical equation can be obtained as above
80 heat flux 50
for the heat transport in a fin. C is the
75 0 dimensionless concentration of the reaction
0 0,02 0,04 0,06 0,08 0,1 component, X the dimensionless coordinate and
length of fin x [m] φL is the Thiele modulus for the plate and a first
order reaction.
Figure 1. Temperature and heat flux depending on the
length of fin. C ′′ − φ L2 ⋅ C = 0
c cosh(φ L X )
A slightly more complicated example C( X ) = =
cK cosh(φ L )
concerns the accompanying heating of a pipe.
The students have to prepare the technical k
drawing with a CAD program. Then this φ L2 = L2 ⋅
Deff
“geometry” can be imported. The next steps
concern the transformation to a solid and the For the estimation of the catalyst effective-
conversion of millimeters to meters. The ness factor there are two possible ways to obtain
temperature field can be calculated for given the calculation of the converted amount (by way
values of heat convection coefficients, heat of analytical solution and similarly by way of the
conduction coefficients, fluid temperatures in numerical solution). The first way is the
both pipes and the ambient temperature. subdomain integration over the reaction rate
Figure 2 shows the temperature field for the (source term).
given situation. The large pipe contains the VK
product and the small pipe the heating fluid.
Minimum and maximum temperatures and the ∫
η= 0
k ⋅ c( x) ⋅ dV K
amount of heat loss can be obtained based on k ⋅ cK ⋅V
these results.
The following second way concerns the 2
calculation of the diffusion flux at the boundary C ′′ + ⋅ C ′ − φS2 ⋅ C = 0
X
of the plate.
n& Diff a system is obtained with two ordinary differen-
η= tial equations.
k ⋅ cK ⋅V
In either case the following formula is C′ = Z
obtained from the analytical solution 2
tanh(φ L ) Z′ = − ⋅ Z + φS2 ⋅ C
η= X
φL For a given value of ∆X there is a difference
between numerical and analytical solution.
A modified Thiele modulus allows for the Figure 4 shows the error in calculation, that
comparison of different geometries. Instead of means the above-mentioned difference between
diameter or height the ratio of volume and outer analytical and numerical solution for the
surface is employed. For plates there is no concentration at X = 1 depending on the step
difference: L = V / O and φ L = φ mod size. The method employed is based on the
2 following formula.
V k
φ mod
2
= ⋅
O Deff C1 , Z1 = f ( X o , Co , Z o )
Figure 3 shows the effectiveness factor
depending on the Thiele modulus for different
5%
shapes of catalyst. The greatest difference
error in calculation
between sphere and plate of about 14% is found 4%
at φ mod = 1.5. This includes the results of 3%
calculations with Comsol Multiphysics for a
2%
cylinder (h = d), a hollow cylinder (h = d = 2 di)
and a monolithic catalyst with rectangular canals. 1%
0%
0 0,02 0,04 0,06
1,0
0,9 plate ∆X
sphere
0,8
cylinder Figure 4. Error in calculation (concentration
0,7
hollow cylinder difference) at X = 1
effectiveness factor
0,6 monolith
0,5
0,4 In the case of exothermic reactions it is
0,3
possible, that the catalyst effectiveness factor is
greater than 1. This is due to the enhancement of
0,2
the chemical reaction by the higher temperature
0,1
inside the catalyst particle. Analytical solutions
0,0
do not exist.
0,1 1 10 100
Figure 5 shows calculated results exemplary
modified Thiele modulus for the decomposition of nitrogen oxide (N2O).
The interaction between the temperature and
Figure 3. Effectiveness factor depending on the Thiele
modulus. concentration fields causes a thin layer with a
high reaction rate. For a lower radius the
temperature is high, but the concentration of the
For the case of a sphere catalyst particle the reactant is nearly zero. And for a greater radius
students have to calculate the concentration the concentration is high, but the temperature is
profile for a first order reaction with the Euler not sufficient for the start of chemical reaction.
method while using Excel. From
In this case there is an effectiveness factor of cooling of work pieces to cooking / freezing of
η = 13.3! That means that the reaction rate is food. Analytical solutions exist only for special
about 13 times faster compared with the case of cases. They consist of an infinite sum of partial
neglecting of all transport resistances (equivalent functions. Simple functions are obtained only for
to a reaction rate under surface conditions). The very short or very long periods of time.
comparison of the calculated data with data The selected example concerns the one-sided
reported in literature indicates good agreement. heating of a wall. The practical background
An indication of the quality of the could be a fired wall of a store-room. How
calculations can be determined by using the rapidly does the temperature of the cold side of
combination of energy and mass balances. This the wall increase? An analytical solution exists
combination results in the Prater number β : for the case of adiabatic behavior of the cold
side. With Comsol Multiphysics the solutions
To − TS (− ∆H ) ⋅ Deff ⋅ cS
β= = can be compared for the adiabatic behavior and
TS λeff ⋅ TS for the case of convective heat transfer to the
The same numerical value is obtained in the cold environment. Figure 6 shows the
above-mentioned example by using the temperature of the cold side of the wall
difference between the surface and center depending on the time for both variants.
temperature as well as by using the properties “Convective” means convective heat transfer to
and the surface concentration. the cold environment (20°C) with a heat transfer
The difficulties in numerical calculation are coefficient of 20 W/m²K.
due to the possible existence of very high In postprocessing the heat flux can be
gradients. That is why we have to start with very obtained at both boundaries - hot and cold -
low surface concentrations. Then the depending on time. Figure 7 shows these data.
concentration can be increased step by step. The absorbed heat flux at the start of the heating
Alternatively the parametric calculation method process is very high and the gradient is high as
can be used. well.
85
30
1400 75
isolation
temperature [K], reaction rate [mol/m³s]
temperature [°C]
65
25 convective
1200 55
concentration [mol/m³]
45
1000 20
35
temperature 25
800
reaction rate 15
15
concentration 0 2000 4000 6000 8000 10000
600
10 time [s]
400
Figure 6. Temperature of the cold side of the wall
5
200 depending on time.
0 0
0 1 2 3 4 5
These data can be exported to Excel and then
radius [mm] calculated for the difference between both heat
fluxes (input and output). The integration of this
Figure 5. Temperature, concentration and reaction flux difference over time with the Simpson
rate depending on radius of a spherical catalyst method results in the accumulated energy inside
particle. the wall.
t
∆Q&
Q(t ) = A ⋅ ∫ dt
4.4 Transient Heating of a Wall A
0
“Transient heating” is an important process
in many application fields: from heating /
10000000
1000000 left
right
100000
heat flux [w/m²]
10000
1000
100
10
1
0 2000 4000 6000 8000 10000
time [s]
Figure 7. Heat flux at both boundaries
On the other hand this result (the
accumulated energy) is obtained through
comparison of the initial temperature and mean
temperature at time t. The mean temperature is
obtained by subdomain integration over T(x) for
a given time.
1
T = ⋅ TdV
V ∫
Q (t ) = m ⋅ c ⋅ (T (t ) − To ) = V ⋅ ρ ⋅ c ⋅ (T (t ) − To )
For our data shown in figures 6 and 7 there is a
difference of about 30%. This is due to the
extremely high gradient of heat flux in the first
seconds. This provides a good starting point for a
discussion with students. What influence does
the subdivision in elements have on the accuracy
of results? Which differences can be tolerated?
5. Conclusions
The aim of the course entitled “Transport
processes” is to deepen and consolidate the
students’ knowledge in the fields of heat and
mass transfer, chemical reaction engineering and
fluid dynamics. This should be combined with
education in the field of numerical calculation of
various processes. Comsol Multiphysics is a very
useful instrument to reach this aim. Even the
problems linked to finite element calculations
help students to understand processes as well as
solution methods.
6. Acknowledgement
I wish to thank Dr. Geike Jr. and Dr.
Pocklington for their support in preparing this
paper.