Modeling of Packed Bed Reactors: Hydrogen Production by the Steam

Reforming of Methane and Glycerol

A. G. Dixon*,1, B. MacDonald1, A. Olm1
1
Department of Chemical Engineering, Worcester Polytechnic Institute, Worcester, MA, USA
*Corresponding author: [email protected]

Abstract: COMSOL Multiphysics was used to The standard approach to modeling the complex
model methane steam reforming (MSR) and particle/tube arrangement in a packed bed is to
glycerol steam reforming (GSR) in a fixed bed employ an effective medium approach with one-
reactor of tube-to-particle diameter ratio N = dimensional flow and lumped transport
5.96. The calculations were made using a steady- parameters. Simplification of the flow and
state 1D heterogeneous model, made up of an estimation of the transport quantities is usually
effective medium reactor tube model coupled to done by experiment and empiricism.
a single spherical pellet model for the source The aim of our research is to use COMSOL
term at each point of the tube. The focus of the Multiphysics to model a tubular packed bed
present paper is to show how to include the reactor where the tube is an effective medium,
effects of mole changes due to reaction on the with the reaction rates obtained by solving a
reforming processes in the absence of heat single pellet model at each point [1]. The
effects. The changes in gas density and velocity reforming reactions studied involve changes in
are accounted for by tracking changes in pressure moles, and a rigorous treatment of the effects of
and mean molar mass along the tube. The this on velocity and conversion is presented. In
implementation of this multi-geometry model in this paper we calculate the impact of the mole
COMSOL Multiphysics is discussed, and two changes in the absence of heat effects, so the
different approaches to accounting for the mole simulations reported here are isothermal. We
changes are compared. illustrate the approach using a one-dimensional
domain for the reactor tube coupled to a 2D
Keywords: Reaction engineering, packed bed, domain for the single particle model [2] in
hydrogen production, methane steam reforming, COMSOL
glycerol steam reforming
2. Model Equations
1. Introduction
The model equations are given here in
The production of hydrogen is important in dimensional form. Axial dispersion is included,
the chemical industry; applications include as is axial variation of pressure along the tube.
hydrotreating and energy conversion by fuel The catalyst pellet was taken to have spherical
cells. The conventional route is by the symmetry.
endothermic steam reforming of methane (MSR)
CH4 + H2O ↔ 3H2 + CO 2.1 Ergun equation for pressure drop
CO + H2O ↔ CO2 + H2
CH4 + 2H2O ↔ 4H2 + CO2 The pressure drop in the packed column is
in a multitubular packed bed, at > 20 bar and described by the differential form of the Ergun
700-800 °C and using high flow rates. With the equation:
increasing use of biodiesel as a renewable fuel,
interest has grown in producing hydrogen by
steam reforming of the excess glycerol (GSR)
C3H8O3 + 3H2O → 7H2 + 3CO2 (1)
which is produced as a side product. This 1 1
endothermic reaction also takes place at high 150 1.75
temperature but at lower pressures (see Table 1
in the Appendix for values).
where is pressure (Pa), is the reactor axial
Packed beds are important in the chemical
coordinate (m), is the average bed voidage,
industries in separations and as catalytic reactors.
is the mass flux (kg/m2s), is the catalyst pellet

Excerpt from the Proceedings of the 2014 COMSOL Conference in Boston

diameter (m), is the fluid viscosity (kg/m·s), After substituting for the first term in equation
and is the fluid density (kg/m3). (3), ci is the only dependent variable; the
equation is solved for i = 1, …, Ns. Then the
2.2 Catalyst pellet species balances equations
, (6)
Concentration profiles were computed for
five components in the MSR case and four
components for GSR. The profiles were based on are used to account for the change in velocity u
the following material balances accounting for due to mole changes and pressure drop.
reaction, convection and Fickian diffusion in the The above approximation provides a
spherical catalyst domain: convenient way to deal with mole changes;
however, the drawback is that the total
∗
concentration is not constrained to satisfy the
∗ ∗ ∗
∑ (2) ideal gas law (or other equation of state). A
different approach [4] follows from writing
where ∗ represents the radial position within the
pellet (m), represents the diffusivity of 1
component in the pellet (m2/s), represents
the concentration of component in the pellet (7)
(mol/m3), represents the stoichiometric
coefficient of component in reaction , and
represents the rate of reaction (mol/m3s). where

2.3 Gas phase species balances , (8)

The concentration profiles in the reactor bulk are the mean molar mass and its derivative, and
gas phase are described by G is the constant mass flux (kg/m2·s). If we write
Fick’s law correctly as
1 (9)
(3)
then for the gas phase species balance we get

where is the cross-sectional area of the

reactor tube (m2), is the molar flow rate of (10)
component (mol/s), is the ratio of the
pellet’s surface area to volume (m-1), is the 1
particle-to-fluid mass transfer coefficient (m/s), for i = 1, 2, …, Ns-1. This is supplemented by
is the concentration of component at the
pellet surface (mol/m3), is the concentration of
component in the gas (mol/m3), and is the 1, , (11)
reactor axial dispersion coefficient (m2/s).
Equation (3) is not in the form required by
COMSOL, as there are two dependent variables 2.4 Boundary conditions
Fi and ci which are related by
At the tube inlet we take
(4)
, (12)
Two approaches can be used to address this
problem. In the first approach [3], we can write, while at the tube exit the outflow condition
assuming constant Ac
0 (13)
1
(5)

Excerpt from the Proceedings of the 2014 COMSOL Conference in Boston

was imposed. For the pellet equations there was
symmetry at ∗ , while at the pellet-gas 2 1.1 . . (17)
interface equating the fluxes gave
The tube model equation requires the species
fluxes from the pellet surface as a source term,
∗
| ∗ (14) and these were made available using coupling
variables [3]. Likewise the pellet equations
require the tube bulk concentrations as boundary
3. Use of COMSOL Multiphysics conditons, and again coupling variables were
used. Figure 2 shows schematically how this was
Our approach follows that of the Packed Bed set up in COMSOL Multiphysics.
Reactor example in the COMSOL User Guide
[3], but with a modified treatment of mole
changes. Two model domains were used, a line
for the 1-D tube and a square for the 2-D (a)
spherical particle model, connected by coupling
variables (Figure 1). The coordinates r* and z
were scaled by x = z/L and y = r*/Rp. In the gas
phase the concentration variable was defined by

∗ (15)

Particle

y
Coupling
Cpi(y,x)
(b)
variables

Tube ci*(x)
0 x 1

Figure 1. Schematic of tube-particle models

Tube pressure drop was calculated by the

Ergun equation using the Coefficient PDE
feature of COMSOL. The equations for diffusion
and reaction in the pellet and convective
dispersion in the tube were handled by use of Figure 2. Coupling scheme to supply tube
COMSOL Multiphysics and the Chemical concentrations to particle fluxes (a) and to supply
Reaction Engineering Module. particle surface fluxes to tube balance (b)
Literature kinetics were available for both the
MSR [5] and GSR [6, 7] reactions. The kinetic The domains were meshed using 100
constants are available in the original references. uniform elements for the 1D tube model, while
Standard correlations were used for the axial in the 2D pellet model the same 100 elements
dispersion coefficient Dea [8]: were used in the x-direction, and 60 elements
were used in the y-direction. The y-direction
0.73 0.5 mesh was distributed so that it was refined near
9.7 (16) the pellet surface, using a geometric sequence
1
with an element ratio of 0.1. A typical MRS
and the particle-gas mass transfer coefficient kg solution used 34,010 degrees of freedom and
[9]: GSR used 25,249.

Excerpt from the Proceedings of the 2014 COMSOL Conference in Boston

4. Results and Discussion

The results first show the consequences of

using the concentration-based approach, then
results for MSR and GSR are presented using the
molar mass method advocated here.

4.1 Approach using equations (5) and (6)

A major drawback of the approach using ci as

dependent variables as in equations (5) and (6) is
that the total concentration is not forced to
satisfy the equation of state. This is illustrated
for MSR in Figure 3.

Figure 4. CH4 concentration inside the catalyst pellet

(y-coordinate is r*/Rp) and along the tube (x-
coordinate is z/L) for MSR

To illustrate the effect of pressure, the MSR

simulation was run twice, once with the pressure
drop computed by the Ergun equation, and once
with pressure constant at the inlet value. A
comparison for CH4 conversion is shown in
Figure 5.

Figure 3. Comparison of ctot for MSR using the two

methods described in this paper

The figure shows that the expression P/RT

decreases along the reactor tube due to the
pressure drop in this isothermal simulation. The
quantity ctot which is calculated as in equation (6)
is seen to continually increase, which is not the
expected behavior on physical grounds.

4.2 MSR simulation

Figure 4 shows the concentration of methane Figure 5. Axial profile of CH4 conversion for MSR
in the two-dimensional pellet domain. Along the comparing actual to zero pressure drop cases
x-coordinate at y = 1 we see that the surface
concentration of CH4 decreases monotonically For both cases, the CH4 conversion increases
down the reactor tube as the reactant is used up. at first due to the initial consumption of reactant
Along the y-coordinate at low values of x we see in the tube. For constant pressure, it then levels
a sharp drop near the pellet surface, due to the out as the reactions come to equilibrium. For the
strong diffusion limitations for this reaction. case with pressure drop, equilibrium conversion
Inside the pellet the reaction reaches equilibrium slowly increases following Le Chatelier’s
and the CH4 concentration is constant. principle due to the lower pressure.

Excerpt from the Proceedings of the 2014 COMSOL Conference in Boston

 more strongly down the reactor tube as the GSR
reaction rate is higher than that of MSR. Along
the y-coordinate we see a very sharp drop near
the pellet surface, due to the even stronger
diffusion limitations in GSR than in MSR. Inside
the pellet the C3H8O3 concentration drops to zero
as the reaction is modeled as being irreversible
and thus complete conversion is possible.

Figure 6. Velocity profile along tube for MSR,

comparing actual to zero pressure drop cases

Figure 6 shows how the linear velocity

changes for the same two cases with and without
pressure drop. The conversion in the MSR case
is relatively low, so in the absence of pressure
effects the velocity change due to the increase in
moles is small, only about 1%. Similarly the
change in mean molar mass is small. The larger Figure 8. Axial profile of C3H8O3 conversion for GSR
change seen for the pressure drop case is about comparing actual to zero pressure drop cases
8-9%, mostly caused by the change in pressure
through ctot in equation (11). For the GSR reaction the pressure drop is
lower than for MSR. Both glycerol conversion
4.3 GSR simulation and mean molar mass are virtually unaffected by
the pressure drop along the tube (Figures 8 and
Figure 7 shows the concentration of glycerol 9). The glycerol conversion reaches 25% in only
in the two-dimensional pellet domain. a 1 m long reactor. The mean molar mass
decreases by 12% in this case.

Figure 7. C3H8O3 concentration inside the catalyst

pellet (y-coordinate) and along the tube (x-coordinate)
for GSR.
Figure 9. Axial profile of mean molar mass for GSR
Along the x-coordinate at y = 1 we see that comparing actual to zero pressure drop cases
the surface concentration of C3H8O3 decreases

Excerpt from the Proceedings of the 2014 COMSOL Conference in Boston

 4. Cropley, J. B.; Burgess, L. M.; Loke, R. A.,
The optimal design of a reactor for the
hydrogenation of butyraldehyde to butanol, ACS
Symp. Ser., 237, 255-270 (1984)
5. Hou, K.; Hughes, R., The kinetics of methane
steam reforming over a Ni/α-Al2O3 catalyst,
Chem. Eng. J., 82, 311-328 (2001)
6. Cheng, C. K.; Foo, S. Y.; Adesina, A. A.,
Glycerol steam reforming over bimetallic Co-
Ni/Al2O3, Ind. Eng. Chem. Res., 49, 10804-
10817 (2010)
7. Iliuta, I.; Iliuta, M.; Radfarnia, H., Hydrogen
production by sorption-enhanced steam glycerol
reforming: sorption kinetics and reactor
simulation, AIChE J., 59, 2105-2118 (2013)
Figure 10. Axial profile of linear velocity for GSR 8. Edwards, M. F.; Richardson, J. F., Gas
comparing actual to zero pressure drop cases dispersion in packed beds, Chem. Eng. Sci., 23,
109-123 (1968)
Changes in linear velocity in these isothermal 9. Wakao, N.; Funazkri, T., Effect of fluid
simulations are due to both changes in pressure dispersion coefficients on particle-to-fluid mass
and to changes in moles due to the stoichiometry transfer coefficients in packed beds, Chem. Eng.
of the reactions. For MSR the changes seen in Sci., 33, 1375-1384 (1978) 
linear velocity were relatively small and were
due almost entirely to the pressure drop. The 7. Acknowledgements
changes in linear velocity shown in Figure 10 for
the GSR reaction, on the other hand, were
Some of the preliminary calculations for
mainly due to the change in moles in the GSR
glycerol steam reforming were made by John E.
reaction. This also explains the similar behavior
Kent in his senior year qualifying project.
in the pressure drop and constant-P cases.
8. Appendix
5. Conclusions
Table 1: Constants used in MSR and GSR models
The effects of mole increases on the gas
velocity and conversion for MSR and GSR were Quantity
simulated rigorously using the variable mean MSR GSR
(f : feed)
molar mass. Results are reported for isothermal
1-D calculations; future papers will present L 10 m 1m
results for 2-D non-isothermal cases.
dp 0.0254 m 0.0200 m
6. References R 0.075692 m 0.0596 m

1. Allain, F.; Dixon, A. G., Modeling of uf 2.62 m/s 2.00 m/s

transport and reaction in a catalytic bed using a kg 0.1378 m/s 0.618 m/s
catalyst particle model, Proceedings of the 2
COMSOL Conference, Boston (2010) Dea 0.033 m /s 0.0663 m2/s
2. Petera, J.; Nowicki, L.; Ledakowicz, S., New ρf 3.87 kg/m3 0.75 kg/m3
numerical algorithm for solving
multidimensional heterogeneous model of the μ 3.78×10-5 kg/m·s 2.74×10-5 kg/m·s
fixed bed reactor, Chem. Eng. J. 214, 237-246
T 1019.05 K 823 K
(2013)
3. COMSOL Multiphysics Model Library, Pf 21.1 bar 2.02 bar
Chemical Reaction Engineering Module, Packed
bed reactor, Version 4.3b, COMSOL AB (2013)

Excerpt from the Proceedings of the 2014 COMSOL Conference in Boston