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.
Email: cdb36@cam.ac.uk
Keywords: Lattice Boltzmann, gassolid reaction, temperature
Abstract
The lattice Boltzmann (LB) method has been applied to model the gassolid 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 gassolid reactions at both the particle and reactor level.
Introduction
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 NavierStokes 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
biomass:
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
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 NavierStokes 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 nonisothermal
flows. Below, the methods used to model the hydrody
namic, mass transfer and reaction processes will be de
scribed.
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)
jLBixl) = 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):
1LB
3 sis
The current LB model recovers the incompressible
NavierStokes 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:
Ax"
v = VLB
at
1\ Ax"
2/ at
1
Ts
3
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 multicomponent,
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)
(12)
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
Here,
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
required.
Experimental
For comparison to LB simulations, particles of iron ox
ide were prepared by spraying water, purified by reverse
osmosis, on to Fe20s powder (SigmaAldrich; 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 (KleerMount, MetPrep)
and crosssectioned 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 012 vol. %b (ABB) and [CO,] in the
range 020 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, USBTemp). 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 fis
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 wellmixed batch reactor. For example,
after streaming for the first order reaction
clPs,LB(X, t)
,lfLE:
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)] = )
(19)
x [f;,,(x, t) (Ps,LB + Es Ps LB)Wi (1 + 4ci uLB)I
Ps,LB Ps,LB
+ + 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.
Results
To verify the hydrodynamics in the D3Q19 lattice Boltz
mann code, Poiseuille flow was simulated. Half
way bounceback 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 Mimes. 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
S0
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 NavierStokes 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:
(k
L~y,LB 2
(20)
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)
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
follows:
Fl0.1
0.01
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
wellmixed. 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
kilol
10.3
113
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
Pre2o,
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
<@co
Pre2os
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 < Lz,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
Prneriment
lattice Boltzmann
01.0
1.0
In1 a.0 O
LU~ 4
LxLB
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 10s 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 24;
Figure 3: Experimental crosssection (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
elled.
Figure 3 shows the experimental crosssection 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 cutoff 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 106 s/time step .
Physical
Value
0.0102
0.0102
0.08
0.0002
0.33
0.61
0.000122
0.6
10320
1023
2.4 x 10' exp(75000/RT)
1.5 x 104
2.5 x 10s
1.6 x 104
1.0 x 104
Lattice
Value
102
102
800
20
0.00833
1/3 1/9 1/2)
17200
0.00485
1/4(7ioo 1/2)
1/4(7ioo 1/2)
1/4(7T 1/2)
1/4(77 1/2)
Symbol
Lx
L,
Lz
d,
p
uo
v
Pco
pFe203
T
kg
Dext
Dint
aet
aint
Units
m
m
m
m
kg/m
m/S
m2/S
mol/m
mol/m
K
s
m2/s
m2/s
m2/S
m2/S
Symbol
Lx,LB
Ly,LB
Lz,LB
dp,LB
pLB
UO,LB
VLB
pCO,LB
pFe203,LB
TLB
ki,LB
Dext,LB
Dint,LB
Olext,LB
Olint,LB
Relaxation Parameter
1/S9 = 0.55
Scaling
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 24. 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 24
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.
550
500
450
400
350
300'
250
200'
to 0 10
Y 20 '20 60
comresionerror
'LB
exp 2
exp 3
exp 4
0 20 40 60 8(
distance, (mm)
LB
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) 
S0.6
0.4
S0.
a 1()4(
exp 2
exp 3
exp 4
500
time, (s)
1000
(b) loso
1040
() 2()( 4()( 6()()
time, (s)
8()( 1()(()
1020
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
maxunum.
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 gassolid reactions has been demonstrated and appears
promising in the modelling of the fundamental transport
processes involved.
Acknowledgements
Financial support from the Engineering and Physical Sci
ences Research Council (EP/FO27435/1) and the Gates
Cambridge Trust is gratefully acknowledged.
References
Chen S, Doolen GD. The lattice Boltzmann method for
fluid flows. Anna. Rev. Fluid Mech. 1998;30:329364.
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:76237630.
Barin I, Knacke O. Therinodynainical properties of inor
ganic substances. Berlin: SpringerVerlag, 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:125.
D'Humidres D, Ginzburg I, Krafezyk M, Lallemand P,
Luo LS. Multiplerelaxationtime lattice Boltzmann mod
els in three dimensions. Phil. Trans. R. Soc. Lond. A.
2002; 360:437451.
Sullivan SP, Sani FM, Johns ML, Gladden LF. Simulation
of packed bed reactors using lattice Boltzmann methods.
Chein. Eng. Sci. 2005; 60:34053418.
Flekkey, EG. Lattice BhatnagarGrossKrook models for
miscible fluids. Phys. Rev. E. 1993; 47:42474257.
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
612.
Discussion
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 nonuniformity 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 (106) s per time step and simulation of
the entire reduction reaction would require t 700 s.
Conclusion
The lattice Boltzmann method has been applied to model
the gassolid 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.
1997;87:115136.
Thiele EW. Relation between catalytic activity and size of
particle. Ind. Eng. Chem. Res. 1939;31:916920.
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:10161029.
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: McGrawHill, 1980.
Zou Q, He X. On pressure and velocity boundary condi
tions for the lattice Boltzmann BGK model. Phys. Fluids.
1997;9:15911598.
Zhavoronkov NM, Aerov ME, Umnik NN. Hydraulic re
sistance and packing density of a layer of grains. Zhur
nal Fizicheskos Khimii. 1949;23:342360 (original text in
Russian).
Freund H, Zeiser T, Huber F, Klemm E, Brenner G, Durst
F, Emig G. Numerical simulations of single phase reacting
flows in randomly packed fixedbed reactors and experi
mental validation. Chem. Eng. Sci. 2003;58:903910.
