Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 14.6.3 - Validation of a lattice Boltzmann simulation of metal oxide reduction in a packed bed using experimental measurements
Full Citation
Permanent Link:
 Material Information
Title: 14.6.3 - Validation of a lattice Boltzmann simulation of metal oxide reduction in a packed bed using experimental measurements Reactive Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Bohn, C.D.
Müller, C.R.
Scott, S.A.
Dennis, J.S.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: Lattice Boltzmann
gas-solid reaction
Abstract: The lattice Boltzmann (LB) method has been applied to model the gas-solid reduction reaction of Fe2O3 to Fe3O4 in a mixture of CO+CO2+N2. A multiple relaxation time LB method is used to model the fluid flow; a single relaxation time method which models the reactant species and temperature as passive scalars was applied to incorporate mass and heat transfer. The reaction of (i) a single sphere and (ii) a packed bed of Fe2O3 was investigated and direct comparisons with experimental measurements are presented. The LB approach satisfactorily captures the shrinking core between Fe2O3 and Fe3O4 in a single particle, the concentration profile along a packed bed and the maximum temperature rise in the packed bed during reduction. A discrepancy in the shape of the temperature front between the model and experiments is observed and is a result of the passive scalar method applied to the temperature. The LB method shows promise as an approach to model gas-solid reactions at both the particle and reactor level.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00361
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1463-Bohn-ICMF2010.pdf

Full Text

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

Validation of a lattice Boltzmann simulation of the reduction of a metal oxide with
experimental measurements

C.D. Bohn*, C.R. Millleri, S.A. Scottt and J.S. Dennis*

Department of Chemical Engineering and Biotechnology, University of Cambridge,
Pembroke Street, Cambridge CB2 3RA, U.K.

t Department of Mechanical and Process Engineering, Institute of Energy Technology
ETH Zurich, 8092 Zurich, Switzerland.

i:Department of Engineering, University of Cambridge,

Trumpington Street, Cambridge CB2 1PZ, U.K.

Keywords: Lattice Boltzmann, gas-solid reaction, temperature


The lattice Boltzmann (LB) method has been applied to model the gas-solid reduction reaction of FeOs to FesO4 in
a mixture of CO+CO,+No. A multiple relaxation time LB method is used to model the fluid flow; a single relaxation
time method which models the reactant species and temperature as passive scalars was applied to incorporate
mass and heat transfer. The reaction of (i) a single sphere and (ii) a packed bed of Fe20s was investigated and
direct comparisons with experimental measurements are presented. The LB approach satisfactorily captures the
shrinking core between FeOs and FesO4 in a single particle, the concentration profile along a packed bed and the
maximum temperature rise in the packed bed during reduction. A discrepancy in the shape of the temperature front
between the model and experiments is observed and is a result of the passive scalar method applied to the temper-
ature. The LB method shows promise as an approach to model gas-solid reactions at both the particle and reactor level.


Packed beds are widely used in industry for separation,
drying, heat exchange and heterogeneous catalytic pro-
cesses. Modelling fluid flow and heat and mass transport
in packed solids, however, is difficult given the complex-
ity of the solid boundaries.
The lattice Boltzmann (LB) method (Chen and No-
ble 1998) provides an alternative to traditional finite vol-
ume or finite difference approaches. The LB method re-
covers the Navier-Stokes equation and is capable of han-
dling complex boundary geometries for incompressible or
slightly compressible flows. In many chemical processes
involving low Re flow and moderate enthalpies of reac-
tion, and therefore, mild temperature excursions, this as-
sumption of incompressible flow is satisfactory.
The need for sustainable production of environmentally
benign fuels such as He is now widely accepted. One
method for the production of pure He for use in low tem-
perature polymeric membrane fuel cells, involves the re-

dox reactions of iron and its oxides (Bohn et al., 2008).
Here, iron oxide is first reduced in a mixture of CO and
He (syngas) produced e.g. by the gasification of coal or

FeOs + 3CO 2 Fe + 3CO,.

Condensation of any water yields a pure stream of CO,
suitable for sequestration. Subsequent reoxidation of the
iron material in steam generates pure Hg:

3Fe + 4HO FesO4 + 4H,.

Complete oxidation of FesO4 to FeOs in steam
is not possible owing to thermodynamic constraints
(Knacke, 1973), so air must be used:

2 FesO4 + 1/2 O, = 3 FegOs.

The species of iron involved are Fe, FeO, FesO4 and
Fe20s. The first reduction reaction involving CO is given

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

are the relaxation rates for the D3Q19 model as in
D'Humi~res et al. (2002) and the distribution f is contin-
uous and describes the number of particles at a particular
node of the lattice, x, at time I with velocity, ci. The
conserved moments are the mass density, pLB, and the
momentum, jLB given by


3 Feg03 + CO 4 2 FesO4 + CO, (4)

where the enthalpy for the reaction as written is
AHYO23 K -44 kJ/mol.
The objective of this study is to simulate this reduc-
tion reaction (4) in an individual particle and in a packed
bed and perform a quantitative comparison to experimen-
tal data.

1 Theory of the lattice Boltzmann method

The lattice Boltzmann (LB) method has several advanta-
geous for modelling reactive fluid flow through porous
media: (i) all computations are local meaning it can be
implemented on several processors in parallel with lit-
tle information exchange at each time step (ii) it useS
linear advection and collision operators, but still recov-
ers the incompressible Navier-Stokes equations in the low
Mach number limit, thereby avoiding the computationally
difficult solution of either the Boltzmann or the Navier-
Stokes equation (iii) advection is carried out in a discrete
velocity space with arithmetic operations to avoid round-
ing errors and (iv) it easily accounts for irregular bound-
aries. The disadvantages of the LB method include the
computational requirements involved and the difficulty
of modelling temperature gradients for non-isothermal
flows. Below, the methods used to model the hydrody-
namic, mass transfer and reaction processes will be de-
In the present work, a 3D multiple relaxation time
(MRT) (Lallemand and Luo 2003) lattice Boltzmann ap-
proach with 19 discrete lattice velocities (D'Humieres et
al., 2002) was used to simulate the hydrodynamic flow of
the gas mixture. In the lattice Boltzmann method, particle
transport is simulated by moving continuous particle dis-
tribution functions, f,, from one lattice node, x, to neigh-
boring nodes, x + ciAt x + ax; particle collisions
are modeled by a linear relaxation to a local equilibrium
distribution, f,"g. Briefly, the lattice Boltzmann equation
with MRT collision operator for modeling hydrodynamics

f (x + ax, t + at) f (x, t) (5)
-M- S[m(x, t) meqX,t)l,

where M is a linear mapping from space F R "W to mo-
ment space MI R W such that

m = Mf, f = M m. (6)

PLB (X, t) = fi (X, t)

jLBix-l) = pOLBLB(l~X~f) = CCili(X.1!

where PO,LB can be considered a measure of the mean
density and is set to unity and ci is the discrete velocity
set. For the present D3Q19 model, the kinematic viscosity
is given by (D'Humieres et al., 2002):


3 sis

The current LB model recovers the incompressible
Navier-Stokes equation up to order O(u B). The rela-
tionship between the kinematic viscosity in lattice units,
VLB, and the kinematic viscosity in physical units, v, is
given by:

v = VLB

1\ Ax"
2/ at


Since (i) CO and CO, were diluted in N2 and (ii) CO and
N, have identical molecular masses of 28 kg/kmol, the
fluid properties, e.g. v and p, were assumed to be equal
to that of the primary component of the mixture, Ng, at
1023 K.
To incorporate mass and heat transfer into the lattice
Boltzmann model, the method of Sullivan et al. (2005)
based on that of Flekkey (1993) for multi-component,
miscible passive scalars on a D3Q7 lattice was used. The
lattice Boltzmann equation for mass or heat transfer was
given by:

f, (x+ Ax, t+ At)- fs(x, t) = [f, (x, t)- feq~x X, 1)
where 7, is the relaxation time for species s, either the
reactant, CO, or the temperature, T Additionally,

fI~ = WiPs,LB [1 + 4(ci uLB)I

where ci are the lattice velocities, the subscript i denotes
the lattice direction, ps,LB is the concentration of compo-
nent s given by


Ps,LB(X,t) =) fi,s,(X, t)

S diag(so, sl,..., s16

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

which recovers the correct APs,LB at each time step as


For comparison to LB simulations, particles of iron ox-
ide were prepared by spraying water, purified by reverse
osmosis, on to Fe20s powder (Sigma-Aldrich; purity >
99.9 wt. %b) while mechanically mixing. The resulting ag-
glomerates were heated at 1173 K for 3 h and then sieved
to cl, +1700, 2046 pm. The BET surface area of the
sintered Fe20s particles was ~ 1.0 m /g (Micromeritics,
Tristar 3000); the porosity, measured by Hg porosimetry
(Micromeritics, Autopore IV) was~ 0.60. Particles
were fixed in an acrylic resin (Kleer-Mount, MetPrep)
and cross-sectioned prior to analysis with an optical mi-
croscope (Leica DM/LM, dark field) at 5 x ml.lani lk~.Il il'l
Image analysis was performed with conventional software
(analySIS Image Processing).
Gas was supplied to the reactor from cylinders of (i)
pure CO,, (ii) 10 vol. %b CO, balance N, and (iii) pure
Ng; (iv) laboratory air was also supplied. The flowrates
of streams (iii)-(iv) were measured using mass flow me-
ters (AWM5 101VN, Honeywell), but with streams (i)-(ii)
rotameters were used. Automatic switching between inlet
gases was performed with miniature solenoid valves. The
composition of the effluent gas was continuously mea-
sured using nondispersive infrared analysers measuring
[CO] in the range 0-12 vol. %b (ABB) and [CO,] in the
range 0-20 vol. %b (ABB).
Two reactors were used in experiments. First, a flu-
idised bed made of recrystallised Alg0s (i.d. 20 mm)
with a perforated plate distributor as described in Bohn et
al. (2010), was used to investigate the reduction of sin-
gle particles of Fe20s. The distributor had 5 holes, each
1 mm dia., four aligned in a square array of side 8 mm,
with one central hole. The fluidised material consisted
of 15 g of AlOs (purity > 99.9 wt. %b) sieved to
<1, +300, -425 pm. Second, a packed bed made
from 316 stainless steel tubing (i.d. 10 mm; total length
450 mm) was used. A stainless steel plate (1 mm thick)
Supported the bed: the plate contained 6 evenly spaced
holes in a circular pitch, each with a diameter of 1 mm.
Stainless steel thermocouples (o.d. 1.5 mm, type K) were
positioned every 30 mm along the packed bed by weld-
ing fittings (Swagelok) to the steel tube. The tempera-
ture of the thermocouples was recorded at a frequency of
1 Hz (Measurement Computing, USB-Temp). The reac-
tors were placed in a tubular furnace to heat to the desired
temperature of 1023 K.
The experimental conditions including the flowrates,
gaS COmposition and reactor used are now given. In all
CaSes, the reduction of Fe20s to FesO4 was studied at
1023 K and 105 Pa, using mixtures of CO+CO,+N, (re-

and w = (1/4, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8). Here, uLB
is the macroscopic velocity from the hydrodynamic lat-
tice Boltzmann model given in Eq. (9). As shown in Sul-
livan (2006) this choice of f~ recovers the advection-
diffusion equation up to O(u B) with the diffusivity of
the passive scalar given by:
1 1
DLB 4 (). (15)

The relationship between the diffusivity in lattice units,
DLB, and the molecular diffusivity in physical units, D,
is given by:

D = DLB 4"

1\ Az?
2/ }at

Diffusion within the particle matrix was set to a value
lower than that in the interstices by adjusting the relax-
ation time, 7,, for the species distribution functions de-
pending on whether a lattice node was solid or fluid to
give Dint < Dxt. A second distribution function fi-s
was used for the temperature, T, with the thermal diffu-
sivity given by:

iYLB 4 (Tr --). (17)
with scaling to physical units as in Eq. (16). Again, the
relaxation time, 7T, was independently adjusted depend-
ing on whether a lattice node was solid or fluid to give
aint < aext. Since (i) temperature differences within
the bed were experimentally observed to be moderate (<
50 K at an operating temperature of 1023 K) and (ii) the
reactant and product species, CO and CO, respectively'
were diluted in N, such that at the inlet [CO]/[CO,]/[N,]
5/15/82 vol. %b the passive scalar approach for both
temperature and concentration was deemed adequate.
To incorporate chemical reactions into the lattice Boltz-
mann model, the approach of Sullivan (2005) is used
which permits product species to be not only diffused, but
also advected at the fluid velocity uLB. Each lattice node
is considered a well-mixed batch reactor. For example,
after streaming for the first order reaction

clPs,LB(X, t)

ki,LBPs,LB(X,t) = Ps,LB(X,t). (18)

Evaluation of f~ in Eq. (13) is then performed using the
density, P,,LB t At Ps,LB t Ts 7P, ,LB in plaCC Of
Ps,LB. This implementation yields from Eq. (12) and (13):

) [f;,, (x + ax, t + at) f;, (x, t)] = )

x [f;,,(x, t) (Ps,LB + Es Ps LB)Wi (1 + 4ci uLB)I
+ + aPs,LB = Ps,LB'

kii,LB C Solid(X).

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

spectively, 5, 15 and 80 vol. %b) during reduction and air
(21 vol. %b Os, 79 vol. %b Ng) during oxidation. From
thermodynamics (Knacke, 1973), these compositions will
yield FesO4 and Fe20s, respectively. For experiments
with the fluidised bed, the total flowrate of gas was fixed at
1.5 x 104 m /s, giving E 3/U,f, 7 for the Alg0s bed
material as determined using the correlation of Wen and
Yu (1966). For experiments with the packed bed the total
flowrate was 1.4 x 10-" m3/s giving Rel, = odz,/[v(1 -
~bed)] 17 where Un 0.61 m/s is the superficial ve-
locity at the inlet at 1023 K and 105 Pa, cl, 2 mm is the
particle diameter and v 1.28 x 104 m /s is the kine-
matic viscosity of the inlet gas which was assumed to be
equal to that for pure N, and calculated from Hirschfelder
(1954) and sbed = 0.41.


To verify the hydrodynamics in the D3Q19 lattice Boltz-
mann code, Poiseuille flow was simulated. Half-
way bounce-back boundary conditions were used at the
walls (He et al. 1997). In the flow direction, periodic
boundary conditions were used. The lattice dimensions
were L~y,LB Loz,LB 1 and 5 < L,2,LB < 53; the kine-
matic viscosity was fixed by choosing sg = is 1.0 as
per Eq. (10). The initial condition waS PLB 1 ULB = ;
and the f; were initialised using f M-imes. The fluid
was accelerated using a uniform body force to approxi-
mate a pressure gradient. The force, Fz,LB, WaS varied to
give a maximum centerline velocity of a2,LB 0.05 in
all cases. The force was incorporated according to Lalle-
mand and Luo (2003). Simulations were run until the con-
vergence criteria of

o 0.5


0 -0.5

0 0.2 0.4 0.6 0.8
normalised velocity, troLB / max(utoLB

Figure 1: Comparison between the analytical solution to
the Navier-Stokes equation and results from the
D3Q19 lattice Boltzmann model for Poisieulle
flow in a channel of width 3, 5, 11, 15, 31 and
51 fluid nodes.

where <1, is the external diameter of the sphere, Dint is
the effective intraparticle diffusivity and k; is the intrinsic
rate constant. Here, a cubic domain with an even number
of nodes, viz. L,2,LB L~y,LB Loz,LB = 2W WaS used.
The sphere was constructed using:


L~y,LB 2


with r {2I.K, 4.K, Y.K, 24.K}, giving spheres with true
lattice diameters of clp,LB { 6, 10, 20, 50, 100}. The
centre of the sphere was not a lattice node. The macro-
scopic fluid velocity was set to uLB(X, t) = 0 Over the
entire domain. A fixed concentration boundary condi-
tion was implemented at the particle surface by specifying
Ps,LB = 1 for all exterior fluid nodes using Eq. (13) after
the streaming step with an equilibrium boundary condi-
tion. Here, Ts = 1 was fixed and therefore Dint was fixed
and ki,LB WaS varied to give p ="" {001 0.1 1, 1, 10}.
The observed rate at steady state was evaluated as:

1 ~r\ n n ~n r\* n~ 1 .1

|< IU,LB(X, t) -- h,~LB(X, 1
z z,~LB(X, t)

1) In
-< 010

was satisfied over the fluid nodes in all cases. Fig. 1 shows
that the simulated values of a, were in good agreement
with the analytical solution. In addition, U,,LB < 1 1,
i: < 101 and PLB = 1.0 f 106 OVer the entire
fluid domain. Similar plots were obtained over a range of
viscosities, 0.51 < 1/.ss < 2.5.
To verify mass transfer and reaction in the D3Q7 lat-
tice Boltzmann model, a first order reaction in a spherical
catalyst pellet was considered. The effectiveness factor
of a catalyst is a ratio of its observed rate of reaction di-
vided by its intrinsic rate of reaction in the absence of
mass transfer limitations. Thiele (1939) showed that for
spherical catalysts the effectiveness factor could be gen-
eralised for a first order reaction to:

ki,LB CPs,LB(X).Solid(X)

prior to the calculation of the equilibrium distribution
where .solicl(x) takes a value of 1 or 0 depending on
whether the lattice node is solid or fluid. The theoretical
rate in the absence of mass transfer limitations was

il = (9 coth(p) 1) p = .

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

in an initial change in the effluent CO by <0.5 vol. %b or
10 %b of the inlet value. After a specified time, the inlet gas
to the bed was switched to inert Ng using solenoid valves
controlled by a timer to prevent further reaction. The bed
was subsequently cooled to room temperature under the
flow of N,; the reacted particles were sieved to separate
them from the inert Alg0s bed material and then cross-
sectioned for analysis with an optical microscope.
For the lattice Boltzmann simulations, a particle with
clp,LB = 50 lattice nodes was initialised using Eq. (22),
giving AJ:r 0.04 mm/lattice node for a particle with
cl = 2 mm. The effective diffusivity of CO was calcu-
lated from Dint = (1/DA + 1/Dxt)-1/7 where Dxt
is the molecular diffusivity of CO in the gas mixture ap-
proximated from Hirschfelder et al. (1980), = = 0.6 is
the porosity of the particle, Dr is the Knudsen diffusivity
approximated from Satterfield (1980), and 7 2.4 is the
tortuosity. A final value of Dint = 2 x 105 m /s was
used, giving at = 1 x 105 s/time step for a relaxation
time of 7, = 1 using Eq. (15). The intrinsic rate con-
stant k; was determined from Bohn et al. (2010) and gave
k; = 2.4 x 107 exp(-75000/RT) 3550 s- or in lat-
tice units ki,LB = 0.0355. Thus, for the current problem
p 20. The total number of moles of FeOs initially
present was determined using the mass density of FeOs
(5150 kg/m3) and porosity of the particles ( = 0.6) as



0.1 1 10
Thiele modulus, #

Figure 2: Effectiveness factor, il, versus Thiele modulus,
9, for a spherical catalyst pellet. The analytical
solution (-) and lattice Boltzmann solutions for
simulations with 6 < <1, < 100 calculated us-
ing a D3Q7 model and a passive scalar method
are shown.

The ratio of the observed rate Eq. (23) to the theoreti-
cal rate Eq. (24)gave the effectiveness factor, il, which
is plotted in Fig. 2 as a function of the Thiele modulus,
p. Noticeably, at p 10 a larger deviation from the
analytical solution is observed for spheres of smaller di-
ameter. The error observed for 9 > 1 can be explained
as follows. First, to increase 9 at a fixed 7, 1 and
dp,LB, ki,LB must be increased as shown by Eqs. (15) and
(21). However, increasing ki,LB Will CVentually result in
O(Ps,LB) = O(ki,LBPs,LB) O(APs,LB). That is the
change in Ps,LB fOT a Single time step will be on the order
of the existing Ps,LB. \iglllit.lIII CrOTS in the approxima-
tion of clPs,LBldti Will TCSult and for extreme cases lead
to negative values of Ps,LB. At p 100, only the sphere
with clp,LB 100 returned positive Ps,LB OVer the entire
domain. Second, as 9 is increased the reaction only oc-
curs near the boundary of the sphere. For large p if the
grid spacing A r/R is not kept less than the width of the
front J~substantial errors as a result of the low
accuracy of the equilibrium boundary conditions will also
be incurred.
Next experimental results and lattice Boltzmann simu-
lations for the reduction of a single, spherical particle of
Fe20s to FesO4 via reaction (1) were compared. For the
experiments, particles of Fe20s with <1, 2 mm were
dropped into the fluidised bed at 1023 K, which was flu-
idised with 1.5 x 104 m /s of gas with a composition
of CO+CO,+N, of 5, 15 and 80 vol. %b, respectively. As
in Bohn et al. (2010), the bed could be assumed to be
well-mixed. In each case, 20 particles corresponding to
~ 0.15 g of Fe20s were introduced to the bed and resulted

5150 kgf 1 kilol FeOs
PoieoIt=U sI x 159.7kgf

(1 0.6) 0.8
x x -
1 1


where the factor of 0.8 comes from the experimentally de-
termined observation that maximum conversions of 80 %b
are achieved (Bohn et al., 2010). The conversion X at
each lattice point was then determined using

X 1

Pre20o t=n

Thus, for X = 1, the true conversion of the particle is
80 %b owing to the factor of 0.8 in Eq. (25). The rate of
reaction as a function of conversion was then


-kipcro (27)

-kipco(1 X)

Since the solid is stationary and not advected, no distribu-
tion function f;,, with streaming and collision steps is re-
quired for the metal oxide and only the macroscopic den-
sity Pre2o, was tracked. As before because of their ease
of implementation, equilibrium boundaries were used af-
ter the streaming step to specify the concentration Pco
0.6 mol/m3 or in lattice unitS PCO,LB 1 in all of the sur-
rounding fluid. The initial condition waS PCO,LB = O fOT

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

simulated using a computational domain of L~r,LB
L~y,LB 102 and Loz,LB 800 lattice nodes. In to-
tal, 297 particles, each with clp,LB 20 lattice nodes
were used. The position of the particles was determined
with a discrete element model (Third et al., 2010) by
dropping spheres into a cylinder and gave a porosity of
~bed 0.41. A schematic of the packed bed comprising
the computational domain is shown in Fig. 4; particles
were positioned between 252 < L-z,LB < 538.
Because the reaction and mass transfer terms do not in-
fluence the fluid flow, hydrodynamic and mass transfer
simulations were conducted separately. First, the steady
state flow field was determined; then the mass transfer
simulation was run independently using the calculated
flow field uLB. For the hydrodynamics, at the inlet a fixed
velocity boundary condition of Zou and He (1997) was
used to give


lattice Boltzmann



In1 a.0 O

LU~ 4


Uo = 0.61 m/s, U0,LB = 0.0083.

No special treatment of corner or edge nodes was neces-
sary. At the outlet, a zero flux boundary condition was
applied after streaming as:

f (., y, Loz,LB, t) f (., Y, Loz,LB 1, t). (29)

The hydrodynamics were verified by comparing the pres-
sure drop AP/L where L is the length of the packed bed
with the correlation of Zhavoronkov et al.(1949) as in Fre-
und et al.(2003) at Re, 17. Agreement in AP/L to
within 4 %b was obtained and confirmed the accuracy of
the flow simulation. A full description of the hydrody-
namic parameters used is provided in Table 1.
Mass transfer simulations were run with the uLB(X, t)
from the hydrodynamic simulations as input. The inlet
boundary used the fixed concentration boundary condition
of Zou and He (1997), giving:

Is = Ps,LB 0f 1 f 2 + 3 + 4 +6) (30)

with Ps,LB 1. At the outlet a zero flux boundary iden-
tical to Eq. (29) was applied. To vary the internal and ex-
temnal molecular and thermal diffusivities, the relaxation
times were varied depending on whether a node was solid
or fluid. Values of the relaxation parameters used for the
mass and heat transfer simulations are given in Table 1.
The reaction was incorporated as per Eq. (19) just prior to
the equilibrium step.
The experimental packed bed consisted of 12 g or ~
1600 particles of Fe20s. The flowrate of the reducing gas
with composition CO+CO,+N, respectively, 5, 15 and 80
vol. %b was set to 1.4 x 10-s m3/s to give Re, 17. The
bed was heated to 1023 K in air and the temperature al-
lowed to stabilise for 10 min prior to reduction. In total
4 cycles of reduction were performed. Figure 5(a) shows
the effluent composition as a function time for cycles 2-4;

Figure 3: Experimental cross-section (left) and lattice
Boltzmann simulation (right) for reduction of
a particle of Fe20s (X 0 ) to FesO4 (X 1)
in CO+CO,+N, (5, 15, 80 vol %b, respectively)
at 1023 K. Scalebar = 500 pm.

solid nodes and PCO,LB 1 for fluid nodes. Simulation
run time was decreased by exploiting the problem's sym-
metry only one eighth of a particle was actually mod-
Figure 3 shows the experimental cross-section and cor-
responding model prediction for the reduction of a single
sphere of Fe20s to FesO4 as a function of time, t. Ini-
tially, at t 0 the particle is entirely Fe20s and X 0.
As the reaction progresses a shrinking core delineating
a front between Fe20s and FesO4 moves progressively
towards the centre of the particle. Qualitatively good
agreement between the model and simulation is observed
by comparing the experimental and model predictions at
t 10, 30 and 40 s. Quantitatively, the thickness of the
FesO4 shell at t 10, 30 and 40 s determined from sim-
ulations using the arbitrary cut-off of X 0.5 gave 280,
480 and 620 pm compared to 260 f 30, 490 f 50 and
590 f 50 pm from experiments. For long times (t oo),
the particle is fully reduced to FesO4 and X = 1.
Next, the reduction of a packed bed of iron oxide was

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

Table 1: Table relating the parameters in physical units to lattice units. Here, Ax: = 0.0001 m/lattice unit, at=
1.37 x 10-6 s/time step .

2.4 x 10' exp(-75000/RT)
1.5 x 10-4
2.5 x 10-s
1.6 x 10-4
1.0 x 10-4


1/3 1/9 1/2)


1/4(7ioo 1/2)
1/4(7ioo 1/2)
1/4(7T 1/2)
1/4(77 1/2)




Relaxation Parameter

1/S9 = 0.55

Lx = Lx,LB Xa
L, = Ly,LB Xa
L, = Lz,LB Xa
d, = dp,LB Xa
P = PLB x 0 33
uo = UO,LB X la
v = VLB X a22 a
pco = PCO,LB x 06
PFe203 PFe203,LB x 06
T = 100 x TLB 1023
kgi = kci,LBla
Dext = Dext,LB X a22 a
Dint = Dint,LB X a22 a
OIext = Olext,LB X a22 a
oint = oint,LB X a22 a

Fig. 5(b) shows the temperature measured by three ther-
mocouples spaced ~ 30 mm apart along the length of
the bed as a function of time for cycles 2-4. The ex-
perimental results demonstrate satisfactory reproducibil-
ity; particularly the breakthrough of CO is nearly identical
over the three cycles studied. In Fig 5, results for cycle 1
were omitted as the reduction of Fe20s to FesO4 in cy-
cle 1 in a previous study (Bohn et al., 2010) was found
to give anomalous results compared to subsequent cycles.
From Fig. 5(a), it is evident that all of the entering CO is
consumed in the reduction reaction for t < 500 s. For
t > 500 s, the CO rapidly increases to its inlet value of
0.6 mol/m The effluent concentration profile thus sug-
gests that reduction in the packed bed is characterized by
a sharp reaction front which travels along the reactor and
ultimately breaks through over 500 < t < 700 s.

The thermocouple measurements also track the
progress of the reaction front. The increase in the reac-
tor temperature is sequentially observed for TI, followed
by T2, then Ts, where TI~, T, and Ts were located progres-
sively further along the length of the bed. The average rise
in temperature for the three thermocouples over cycles 2-4
was AT = 21.0 -E 1.5 K.

To compare the time dependent experiments with
the instantaneous concentration and temperature profiles
from the packed bed, simulations which vary with posi-
tion, z, a change in reference frame based on the speed of
the reaction front was necessary. Using the distance be-
tween thermocouples (~ 30 mm) and the times between
temperature maxima, Atl2, Aats and Atls, the speed of

Figure 4: Schematic diagram of the packed bed consist-
ing of 297 particles. The computational domain
in lattice units was Lx2,LB Ly,LB 102,
Lz,LB 800; for each particle, d, 20 lat-
tice umits.









to 0 10

Y 20 '20 60



exp 2
exp 3
exp 4

0 20 40 60 8(
distance, (mm)

exp 2 T,

0 20 40 60 8C
distance, (mm)

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

(a) -



a 1()4(

exp 2
exp 3
exp 4

time, (s)


(b) loso


() 2()( 4()( 6()()
time, (s)

8()( 1()(()


Figure 5: (a) Concentration of CO in the effluent gas
as a function of time for the reduction of a
packed bed of 16 g of Fe20s to FesOs using
CO+CO,+N, (5, 15, 80 vol %b, respectively) at
1023 K. (b) Corresponding temperature excur-
sion during reduction as a function of time mea-
sured by thermocouples positioned in intervals
of ~ 30 mm along the packed bed.

Figure 6: (a) Concentration profile for the reduction of a
packed bed of Fe2 0s to Fes Og in CO+CO, +N,
(5, 15, 80 vol %b, respectively) at 1023 K ob-
tained by rescaling the abscissa in Fig. 5 by
the speed of the front, <~L/di 0.15 mm/s.
(b) Corresponding temperature along the bed
during reduction as measured by thermocouple,

the reaction front was measured to be:

<1 = 0.154 -E 0.01 unn/s. (31)

To convert the time dependent profiles given in Fig. 5 to
profiles with respect to reactor position, z, the absicissa
in Fig. 5 were multiplied by -dzL/di from Eq. (31) and
offset so as to intersect the model predictions at half the
Figure 6(a) shows the experimental effluent concentra-
tion as a function of position, z, as well as the front sim-
ulated using the lattice Boltzmann method. The value of
concentration as a function of z from the LB simulation
was determined by averaging pco over all nodes, both
solid and fluid, on the inside of the cylinder and convert-
ing to physical units as outlined in Table 1. From Fig. 6(a)
the modelled front is sharper than the experimentally de-
termined front.
Figure 6(b) compares the instantaneous temperature
along the reactor with the results from thermocouple Ti
from Fig. 5. The total temperature rise calculated using
the LB method by averaging T over all nodes interior to

the tube, including both solid and liquid nodes, was 24 K,
in good agreement with the 21 K observed in experiments.
Comparing the temperature fronts, however, a noticable
difference in the temperature profiles is observed at the
end of the reactor. In experiments, the temperature in the
bed beyond the front is approximately equal to that at the
inlet, and the front has a Gaussian shape. In the LB sim-
ulations the temperature increases at the front and then
remains high along the bed until the exit. The reason for
the persistently high temperature in the LB simulations
towards the exit has to do with the method of incorporat-
ing temperature in to the model. Because T is a passive
scalar that is created as the reaction proceeds and by defi-
nition must be conserved, the temperature at any point be-
yond the front must be approximately equal to that at the
front for no accumulation. In other words, no provision
for the destruction of T through, for example, the heating
of solids of a fixed heat capacity is possible. Neverthe-
less, the maximum temperature rise is correctly captured
by the passive scalar method.

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

observed for the concentration profile in the packed bed.
The use of a passive scalar approach to model temper-
ature correctly predicted the maximum temperature rise
in the packed bed, but failed to capture the true shape of
the temperature front. Further development of advanced
thermal lattice Boltzmann methods is therefore necessary.
Overall, the extension of the lattice Boltzmann method
to gas-solid reactions has been demonstrated and appears
promising in the modelling of the fundamental transport
processes involved.


Financial support from the Engineering and Physical Sci-
ences Research Council (EP/FO27435/1) and the Gates
Cambridge Trust is gratefully acknowledged.


Chen S, Doolen GD. The lattice Boltzmann method for
fluid flows. Anna. Rev. Fluid Mech. 1998;30:329-364.

Bohn CD, Milller CM, Cleeton JP, Hayhurst AN, David-
son JF, Scott SA, Dennis JS. Production of very pure hy-
drogen with simultaneous capture of carbon dioxide us-
ing the redox reactions of iron oxides in packed beds. Ind.
Eng. Chein. Res. 2008;47:7623-7630.

Barin I, Knacke O. Therinodynainical properties of inor-
ganic substances. Berlin: Springer-Verlag, 1973.

Lallemand P, Luo LS. Theory of the lattice Boltzmann
method: acoustic and thermal properties in two and three
dimensions.Phys. Rev. E. 2003;68:1-25.

D'Humidres D, Ginzburg I, Krafezyk M, Lallemand P,
Luo LS. Multiple-relaxation-time lattice Boltzmann mod-
els in three dimensions. Phil. Trans. R. Soc. Lond. A.
2002; 360:437-451.

Sullivan SP, Sani FM, Johns ML, Gladden LF. Simulation
of packed bed reactors using lattice Boltzmann methods.
Chein. Eng. Sci. 2005; 60:3405-3418.

Flekkey, EG. Lattice Bhatnagar-Gross-Krook models for
miscible fluids. Phys. Rev. E. 1993; 47:4247-4257.

Sullivan SP. Lattice Boltzmann development for chemical
engineering applications. PhD Thesis. University of Cam-
bridge, Department of Chemical Engineering. 2006.

Wen CY, Yu YH.m A generalized method for predicting
minimum fluidization velocity. AIChE J. 1966;12:610-


In the comparison of the experimental concentration and
temperature profiles with those predicted from simula-
tions, the response time of the infrared analyser and
sampling line for the CO signal and the response time
of the thermcouples and acquisition system for T must
be considered. The response of the analyser and sam-
pling line was estimated by flowing the reaction mix-
ture, CO+CO,+N, with 5, 15 and 80 vol. %b, respec-
tively, through an inert bed of 16 g of Alg0s particles of
d, +1400, -1700 pm heated to 1023 K. A step change
in the inlet concentration was approximated by automati-
cally switching the solenoid valves to a flow of pure No.
By measuring the decay in the CO signal for the flowrate
used in the packed bed experiments, the response time
of the analyser and sampling line was determined to be
~ 1 s. Deconvolution of the signal did not therefore alter
the shape of the front observed in Fig. 5(a) since the front
took ~ 200 s to break through. Similarly, experiments to
approximate the response time of the thermocouples and
acquisition system were performed by rapidly switching
the thermocouple between boiling water and an ice water
bath and vice versa whilst vigorously stirring. The re-
sponse time of the thermocouples was determined to be
~ 0.3 s and did not, therefore, alter the temperature re-
sults. Since deconvolution did not influence temperature
and concentration profiles, in Fig. 5(a) and (b), the mea-
sured signals were shown without deconvolution.
In the simulation, as the front moves down the bed
small deviations from the true steady state profile will be
present owing to the non-uniformity of particle locations
and fluid flow; however, in general the instantaneous tem-
perature and concentration profile allow quantitative com-
parison with experiments. By rescaling the experimen-
tal results using the characteristic front speed, simulation
of the reduction of the entire packed bed can be avoided.
This is particularly important for the current simulations
where at O (10-6) s per time step and simulation of
the entire reduction reaction would require t 700 s.


The lattice Boltzmann method has been applied to model
the gas-solid reduction reaction of Fe20s to FesO4 with
CO+CO,+N2. Comparison between the the LBM model
and experimental results was achieved by considering: (i)
the conversion profile of a single particle of Fe20s, (ii)
the concentration profile for the reduction of a packed bed
of particles of Fe20s and (iii) the temperature profile for
the reduction of a packed bed. For the reduction of a sin-
gle particle, the LB method correctly captured the shrink-
ing core between Fe20s and FesO4 phases. Satisfactory
agreement between the model and experiments was also

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

He X, Zou Q, Luo LS, Dembo M. Analytic solutions
of simple flows and analysis of nonslip boundary condi-
tions for the lattice Boltzmann BGK model. J. Stat. Phys.

Thiele EW. Relation between catalytic activity and size of
particle. Ind. Eng. Chem. Res. 1939;31:916-920.

Bohn CD, Cleeton JP, Milller CM, Davidson JF, Hayhurst
AN, Scott SA, Dennis JS. The kinetics of the reduction of
iron oxide by carbon monoxide mixed with carbon diox-
ide. AIChE J. 2010;56:1016-1029.

Hirschfelder JO, Curtiss CF, Bird RB. Molecular theory of
gases and liquids. New York: John Wiley & Sons, 1954.

Satterfield CN. Heterogeneous Catalysis in Practice. Lon-
don: McGraw-Hill, 1980.

Zou Q, He X. On pressure and velocity boundary condi-
tions for the lattice Boltzmann BGK model. Phys. Fluids.

Zhavoronkov NM, Aerov ME, Umnik NN. Hydraulic re-
sistance and packing density of a layer of grains. Zhur-
nal Fizicheskos Khimii. 1949;23:342-360 (original text in

Freund H, Zeiser T, Huber F, Klemm E, Brenner G, Durst
F, Emig G. Numerical simulations of single phase reacting
flows in randomly packed fixed-bed reactors and experi-
mental validation. Chem. Eng. Sci. 2003;58:903-910.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs