Simulation of non-catalytic gas–solid reactions: application of

grain model for the reduction of cobalt oxide with methane

B. Khoshandam, R. V. Kumar and E. Jamshidi

A mathematical model based on grain model, is carbon dioxide. The results show good agreement for
presented. The model considers dynamic changes of gas experimental conditions where both chemical kinetics
and solid concentrations and temperature in the pellet of and product layer diffusion are significant in the first
the solid reactant. In addition, structural changes of reaction and the reaction obeys a non-linear chemical
solid reactant (porosity and grain size changes) and kinetics in the second reaction.
diffusion of gaseous reactants through the product layer B. Khoshandam (E-mail: bkh@engineer.com) and R. V.
are considered. Effective diffusion coefficient is Kumar (for correspondence – E-mail: rvk10@cam.ac.uk) are
assumed to change with porosity in the pellet as a at the Department of Materials Science & Metallurgy,
quadratic function in the mass balance equation. University of Cambridge, New Museums Site, Pembroke
Effective heat capacity and conductivity are also Street, Cambridge CB2 3QZ, UK; E. Jamshidi (E-mail:
modelled as functions of porosity in the heat balance jamshidim@yahoo.com) is at the Department of Chemical
equation. The finite volume method has been selected Engineering, Amir-Kabir University (Tehran Polytechnic),
for solving mass and heat balance equations in the pellet Hafez Avenue, Tehran, Iran.
and the computations are done by a program developed © 2005 Institute of Materials, Minerals and Mining and
in MATLAB. The results of the program show good Australasian Institute of Mining and Metallurgy. Published
agreement in comparison with analytical solutions for by Maney on behalf of the Institutes. Manuscript received 28
both isothermal and non-isothermal systems. The model January 2004; accepted in final form 23 September 2004.
was used to simulate the reduction of cobalt oxide with Keywords: Simulation, non-catalytic gas–solid reactions,
methane and also the carbon gasification reaction with grain model, reduction of cobalt oxide with methane

INTRODUCTION that may form at the surface of each grain will offer
Non-catalytic gas–solid reactions are commonly increasing resistance to diffusion with time. The grain
encountered in many chemical and metallurgical model has been extensively employed in modelling of
processes. Reduction of metal oxides with reducing non-catalytic gas–solid reactions and applied to a
gases is a typical example in extractive metallurgy that number of systems such as reduction of nickel
is based on non-catalytic gas–solid reactions. A oxide,6–8 hydrofluorination of UO2,9 and reaction of
general type of non-catalytic gas solid reaction may be zinc oxide with hydrogen sulphide.10
represented as: In the present work, the grain model was applied
for the reduction of cobalt oxide pellets with methane.
A(gas) + bB(solid) → cC(gas) + dD(solid) (1)
A program was developed in MATLAB software that
There have been many papers published about is able to solve mass and heat balance equations in
simulation of such reactions in the last three decades. pellets considering both structural changes in the
Different simulation techniques and models have been pellet (such as grain radius and porosity), and
proposed and a good review is presented in the variations in the effective diffusion coefficient. In
literature.1–4 order to reduce computation time, matrix operation
Gas–solid reactions are normally carried out facilities in MATLAB have been used. These
between a gas phase and a collection of solid particles operations can minimise the number of loops applied
contained in a pellet. During gaseous diffusion in the program. The mass and heat balance equations
through a pellet, a chemical reaction may also occur. were discretised using finite volume method and were
The grain model introduced by Szekely et al.5 is solved in the implicit form. Both Patisson et al.11 and
fundamentally based on this idea. The grain model Liu and Bhatia12 used the finite volume method and a
also referred to as the particle-pellet model is based on semi-discrete Petrov-Galerkin finite element method,
the assumption that grains in a pellet are surrounded respectively; however, porosity change, grain product
by pores through which the gas must diffuse to arrive layer resistance, and effective diffusion coefficient
at a sharp interface between the grain particle and the variation as a function of porosity were not
gas phase for reaction to take place. A solid product considered.

C10 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 DOI 10.1179/037195505X39058
 Khoshandam et al. Simulation of non-catalytic gas–solid reactions

MATHEMATICAL MODEL swell. In other words, the grain dimensions are

assumed to change due to the differences in the
Principle molar volume of the products and reactants.
In the present model, we consider a porous pellet (xii) The porosity can change in the pellet. Change in
immersed in a stream of gas mixture containing the porosity is related to that of grain size and
reducing agent. The reducing agent diffuses from the volume of the grain.
bulk of the gas mixture to the porous pellet surface (xiii) Diffusion of reactant gas through the product
and then diffuses into the pellet. During the diffusion layer around the grains is considered.
of gas into the pellet, a chemical reduction reaction is (xiv) Effective diffusivity of gas through the pellet can
deemed to be occurring. The following assumptions change.
and considerations are made in the development of (xv) Combination of the ‘pellet-shape factors’ with
the mathematical model: ‘grain-shape factors’ gives rise to nine possible
(i) The reaction is irreversible and is assumed to be geometrical situations in the most general case.
equimolar or non-equimolar if the gaseous
reactant A were present as a dilute component Equations
in an inert gas. The present model is developed for unsteady state in one-
(ii) The reaction takes place at a constant pressure. dimensional geometry. Independent variables are
(iii) Rate of reaction is defined as a power equation or position in the pellet and time. Other parameters such as
Langmuir-Hinshelwood equation. Many gas–solid concentration of gas through the pellet, temperature,
reactions involve strong adsorption of gas A on the porosity, effective diffusivity, effective special heat
surface of the solid reactant B which can result in capacity, effective heat conductivity, radius of the grain,
fractional order kinetics (Langmuir-Hinshelwood rate of reaction, radius of unreacted core of the grain,
model) or in non-first order power law kinetics concentration of reactant gas on reaction site in the
(Sohn & Szekely model). pellet and finally conversion are dependent variables.
(iv) The regimen can be transient or steady state and
applies to both mass and heat transfer. Conservation equations
(v) The reaction can be endothermic or exothermic. Mass balance equation in pellets – The mass balance
Every chemical reaction is accompanied by heat equation for the reducing gas agent A, based on
and only under special circumstances is this Accumulation = Input – Output + Generation, can be
negligible. expressed as:
(vi) The pellet is porous and has regular geometrical 2 _ fC A i r F -1
g

shape such as sphere, cylinder and flat plate. The = F1 - 1 2 c z F p -1

D e 2C A m + Fg _1 - f 0 i c F f G
2t z p
2z 2z r0 g

balance equations can be derived for different

(2)
pellet shapes using a pellet shape factor, Fp
defined as 1, 2 and 3 for flat plate, cylindrical where all terms are defined in Table 1.
and spherical pellets, respectively.5 The initial and boundary conditions for different
(vii) The pellet is made from grains by compression. pellet geometries (Fig. 1) are:
(viii) The grains are non-porous and have regular
at t = 0, CA = CA0 (initial uniform concentration)
geometrical shapes such as sphere, cylinder and flat
(3)
plate. The source terms for balance equations can be
at z = 0, 2 C A
= 0 (symmetric condition at
derived for different grain shapes using a grain 2z
shape factor, Fg, defined as 1, 2 and 3 for flat plate, the centre of pellet) (4)
cylindrical and spherical grains, respectively.
at z = rp, D e 2C A = k c _ C Ab - C As i
(ix) The grains are so small that temperature 2z
(mass flux equality) (5)
gradients in them can be neglected.
(x) There is diffusional resistance outside the pellet where the terms are defined in Table 2.
from bulk of the gas flow to surface of the pellet.
(xi) During the reaction, the grains can shrink or Mass balance equation in grains – The mass balance

1 Different pellet geometries

Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 C11
Khoshandam et al. Simulation of non-catalytic gas–solid reactions

Table 1 Terms for mass balance equation in pellet Table 2 Terms for initial and boundary conditions of mass
balance equation in pellet
Term Definition
Term Definition
ε Porosity of pellet
CA Concentration of reducing agent A in pellet CA0 Initial concentration of reducing agent A in pellet
t Independent variable, time rp Pellet radius for spherical and cylindrical pellets and
z Independent variable, position in the pellet, radial half-thickness for flat pellets
position in spherical and cylindrical pellets and axial kc Convective mass transfer coefficient
CAb Concentration of A in bulk of gas stream
position in slab-like pellets
CAs Concentration of A on surface of pellet
Fp Pellet-shape factor, 1, 2 and 3 for infinite slab, long
cylinder and sphere, respectively
De Effective diffusion coefficient of A in pellet Table 3 Terms for mass balance equation through the
Fg Grain shape factor, 1, 2 and 3 for infinite slab, long product layer
cylinder and sphere, respectively
ε0 Initial porosity of pellet Term Definition
rc Radial position of reaction front in grain
CAi Concentration of A on reaction front
r0 Initial radius of grain DAp Diffusivity of A in the product layer around grains
fG Rate of surface reaction rg Radius of grain

equation for the solid reactant can be transformed to where ρD, MD, d and εD are density, molecular weight,
an equation for radius of the unreacted core of the stoichiometric coefficient and porosity of solid
grains in terms of Accumulation = Input – Output + product D, respectively.
Generation where Input – Output = 0:
drc =- M B bf Heat balance equation in a pellet – The heat balance
tB G (6)
dt equation for the pellet can be written in terms of
where MB, ρB and b are molecular weight, density and Accumulation = Input – Output + Generation, as
stoichiometric coefficient of solid reactant B, respectively. follows:
The initial condition for this equation is: 2 _ C pe T i
at t = 0, rc = r0 (7) =
2t
1 2 zF r F -1 g

Mass balance equation through the product layer K e 2T m + Fg _1 - f 0 i c F f G DH

c p - 1

z F - 1 2z
p
2z r0 g

around grains – the mass balance equation for the (12)

reducing gas agent through the product layer in terms where the terms are defined in Table 4.
of Accumulation = Input – Output + Generation The initial and boundary conditions are:
(where accumulation and generation terms are zero
and mass flux = reaction rate on surface): at t = 0, T = Tb (initial uniform temperature
For flat grains (grain-shape factor, Fg = 1): in pellet) (13)
fG at z = 0, 2T = 0 (symmetric condition at the
C Ai = D _ rc - rg i + C A (8) 2z
Ap centre of pellet) (14)
For cylindrical grains (Fg = 2):
rc f G at z =rp, K e 2T = h c _ T b - T s i + fv
u _ T w4 - T s4 i
C Ai = D ln c rrgc m + C A (9) 2z
Ap
(heat flux at the surface) (15)
For spherical grains (Fg = 3):
where the terms are defined in Table 5.
rc f G _ rc - rg i
C Ai = D rg + CA (10)
Ap
Auxiliary equations
The terms are defined in Table 3.
The porosity changes in the pellet can be related to
changes in the grain size and the shrinkage–swelling of
Change in grain size – The mass balance equation for
the overall pellet size:
the solid reactant of a grain can be easily simplified to Fg

the following equation for describing the change in rg V p0

f = 1 - _1 - f 0 i d r0 n d V n (16)
p
grain size:
J rF - rF N1 F g

K d _0 g
c i tB M D
g

F O Table 5 Terms for initial and boundary conditions of heat

rg = + rc g
(11)
K b 1 - fD tD M B O balance equation in pellet
L P
Table 4 Terms for heat balance equation in pellet Term Definition

Term Definition Tb Initial uniform temperature

hc Convective heat transfer coefficient
CPe Effective specific heat capacity of pellet Ts External surface temperature of pellet
T Temperature fu Emissivity of the pellet
Ke Effective thermal conductivity of pellet σ Stefan-Boltzmann constant
∆H Heat of reaction Tw Reactor inner wall temperature

C12 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114
 Khoshandam et al. Simulation of non-catalytic gas–solid reactions

where VP0 and VP are initial pellet volume and variable Table 6 Terms for the reaction rate laws
pellet volume, respectively.
Term Definition
Intraparticle diffusion coefficient is defined based
on the random pore model of Wakao and Smith:13 k Reaction rate constant
2

D e = D e0 c ff0 m
CB Concentration of solid reactant B
(17)
n, m Constant values
where De0 is initial effective diffusion coefficient of gas k1 A constant value
reactant in the pellet. k’ Frequency factor
Extent of reaction or conversion about non-porous Ea Activation energy of chemical reaction
grains is defined as local conversion and is the ratio of R Ideal gas constant
mass of solid product D produced at each time to
mass of solid product if the grain converts to solid
nd, the number of solid grains per unit volume of
product completely; so for grains:
Fg porous solid, is calculated from:
x = 1 - c rr0c m (18) 3 _1 - f i
nd = (26)
4rrg3
where x is conversion of the grain and is a function of
local place in the pellet and time. However, conversion The molecular diffusivity is obtained from the
of pellet is defined as overall conversion and is Chapman and Enskog equation.15
calculated from local conversion integration over the The gas viscosity for each component in the bulk of
pellet as: rp the gas can be calculated and the bulk viscosity can be
# z F - 1 xdz
p
obtained using ‘rule of mixtures’.
X = 0 rp (19) The effective specific heat capacity and thermal
# z F - 1 dz p
conductivity of the pellet, respectively, can be
0
where X is overall conversion in pellet and is a obtained from:
function of time and rp is radius of the pellet for
C p = _1 - f i C p + fC p
e s g (27)
spherical and cylindrical pellets and half-thickness for
flat pellets. 2
K e = _1 - f i K s + f 2 K g (28)
Expression of the reaction rate law where CPs and CPg are the specific heat capacity of
The surface reaction rates can be defined as a power solid and gas, respectively, in the pellet and Ks and Kg
equation: are the thermal conductivity of solid and gas in the
pellet, respectively.
f G = kC An C Bm (20)
The convective mass transfer coefficient, kc, can be
and Langmuir-Hinshelwood equation is given by: calculated from the equation used for Sherwood number
(Sh = kc dp/DAM). In general, the dimensionless
f G = kC A (21)
1 + k1 C A Sherwood number can be expressed as a function of
The Arrhenius equation is obtained by considering Reynolds number (Re = ρf u∞ dp/µf) and the Schmidt
temperature effect on the reaction rate as: number (Sc = µf/ρf DAM) with a specific empirical
Ea equation being valid for a given geometrical situation
k = k le - RT (22)
and for a given range of parameters. For example, for a
where the terms are defined in Table 6. laminar flow (0 < Re < 200) past a spherical pellet, the
correlation of Ranz and Marshall,16 is:
Expression of parameters
Sh = 2 + 0.6Re 1 2 Sc 1 3
(29)
The effective diffusivity may be obtained with the aid
of the Bosanquet interpolation formula: In these equations, µf, ρf and u∞ are viscosity, density
-1 and velocity of gas stream, respectively, and dp is
f
De = c x m c D1AM + D1AK m (23) diameter of pellet.

where DAM and DAK are molecular diffusivity and Table 7 Terms for general equation for transport of a scalar
Knudsen diffusivity of gaseous reactant A in the pellet parameter Φ
and τ is tortuosity coefficient of the pellet.
Term Definition
The Knudsen diffusivity is calculated from:
Φ A scalar parameter such as concentration or temperature
1 2
CV Control volume
DK = 4 c r
8RT ko (24)
3 MAm ∆t Time step
dV Volume element
where k0 is calculated from the ‘dusty gas model’ of
dA Surface element
Mason et al.:14
n Outward normal vector to surface element dA
Γ
k 0- 1 = c 128 m n d rg2 c1 + r m
Diffusion coefficient
(25)
9 8 SΦ Source term

Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 C13
Khoshandam et al. Simulation of non-catalytic gas–solid reactions

In an analogous manner, the convective heat temperature as given by Equations (3) and (13). In
transfer coefficient also can be calculated from a each time level, all dependent parameters are also
similar equation used for laminar flow past a sphere: renewed. Next, the solution T and CA are assigned to
T0 and CA0 and the procedure is repeated to progress
Nu = 2 + 0.6Re 1 2 Pr 1 3
(30)
the solution by a further time step.
where Pr is the Prandtl number (Pr = CPf µf/Kf) and
Nu is the Nusselt number (Nu = hc dp/Kf). In these
equations, CPf and Kf are the specific heat capacity and PROGRAM SPECIFICATIONS
thermal conductivity of the gas stream, respectively. The program developed in MATLAB included 12
The experimental conditions are so chosen, in our sections.
work, that Sh and the Nu both tend to infinity, (i) Data input. The data needed for the program
implying rapid mass and heat transfer by convection included the state of the governing equation
in the gas phase. In this manner, the model can be (steady, unsteady isothermal or non-isothermal
tested for selected limiting conditions. state), boundary type (constant value at boundary
[Dirichlet] or a relation between first derivative and
value of parameter at boundary [Neumann]), grain
NUMERICAL METHOD and pellet specifications (weight, geometry, shape
The finite volume integration of partial differential factors), operating conditions (temperature,
Equations (2) and (12), over a control volume (CV) pressure, volume percentage of reducing agent in
with a further integration over a finite time step is the carrier gas, total flow rate, reactor diameter),
applied for the discretisation process of the equations. physical properties of carrier gas, reaction
For transport of a scalar parameter Φ, we obtain:17 components and products (molecular weight,
t + Dt
density, specific heat capacity, heat conductivity),
#d # 2 tU i dt dV = chemical reaction information (stoichiometric
2t _ n
CV t
coefficients), and diffusivity of gaseous reactant
t + Dt t + Dt
through the product layer. In addition for the
# d # n. _ CgradU i dA n dt + # # S U dVdt (31)
t A t CV numerical solution, number of increments in space
where the terms are defined in Table 7. domain, time step and final time should be given
The accumulation term has been discretised using a to the program.
first order backward differencing scheme. Indeed, for (ii) Necessary parameters calculation. In this
the diffusion term, we apply central differencing. The section, parameters needed for running of the
source term, SΦ, can be approximated by means of a program are calculated. These parameters
linear form: consist of: initial porosity, effective heat capacity
and effective conductivity of the pellet,
# S U dV = SU DV = S u + S p U p (32)
CV tortuosity factor, bulk concentration of reducing
where SU is the average of source term in control gas, viscosity of gas mixture, effective diffusivity
volume, Su and Sp are constant values and Φp is of reactant gas through the pellet, Sherwood
parameter Φ in node p inside the control volume. and Nusselt numbers.
The explicit discretisation of equations is very (iii) Chemical reaction rate function. Two rate models
expensive to improve spatial accuracy because can be selected, power equation and Langmuir-
maximum possible time step has to be drastically Hinshelwood equation.
reduced. The two schemes, Crank-Nicolson and fully (iv) Solution of mass balance equation in pellet. In
implicit, have been applied for the discretisation this section, all matrices coefficients according
process. For a general control volume with left-side to the finite volume technique are made and
and right-side face positions of Xw and Xe, then algebraic sets of equations are solved using
respectively, in the Crank-Nicolson scheme, where De matrix operation facilities in MATLAB.
or Ke is constant and a uniform grid spacing, ∆X, we (v) Solution of reaction front equation or mass
have to use time steps as: balance equation for solid reactant. In this
2f _ DX i_ X eF p +1
- X wF + 1i
p section, radius of unreacted core in grain are
Dt < (33) calculated in each time step.
D e _ X eF + X wF i_ Fp + 1i
p p

(vi) Solution of heat balance equation in pellet. In

for Equation (2) and: non-isothermal systems, heat balance equations
2C P _ DX i_ X eF
e
p +1
- X wF + 1i p are solved similar to that used in mass balance
Dt < (34) equations.
K e _ X eF + X wF i_ Fp + 1i
p p

(vii) Calculation of conversion. In the grain and pellet,

for Equation (12). The fully implicit scheme is also the conversion is calculated.
unconditionally stable for any size of time step. (viii) Calculation of grain radius. In this section,
Equation (6) also is discretised using a first order radius of grains are calculated.
backward differencing scheme. After the discretisation (ix) Solution of mass balance equation in grain
process, a system of algebraic equations are solved at product layer. Concentration of gas reactant on
each time level. The time marching procedure starts interface between unreacted core and product
with a given initial fields of concentration and layer around grain are calculated.

C14 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114
 Khoshandam et al. Simulation of non-catalytic gas–solid reactions

a b

2 (a) Effect of grain radius on conversion–time curve. (b) Effect of grain shape factor on conversion–time curve, Fg
= 1 (slab grains), Fg = 2 (cylindrical grains) and Fg = 3 (spherical grains)

(x) Calculation of porosity of pellet. temperature range was 800–950°C. The details of
(xi) Parameter refinement. In this section, all parameters experimental conditions have been described by
such as source term of equation and effective Khoshandam.18 Frequency factor, k′ = 1755·5 m/s and
diffusion coefficient are renewed. The procedure of activation energy, Ea = 155·6 kJ/mol for this reaction
sections (i) to (iv) is repeated for the next time step. were derived from kinetic studies.
(xii) Result visualisation and report file generation.
Finally, in this section, conversion-time plot, Effect of grain size specifications
other necessary 2- or 3-dimensional plots and The effect of grain size has been evaluated by the
report files are generated. model for a slab-like pellet made up of spherical
grains, Fp = 1, Fg = 3 (Fig. 2a). As the radius of the
grain increases, the time needed for completing the
RESULTS OF THE MODEL reaction is increased. For a slab pellet with 3 mm
In this section, effects of changing several important thickness (rp = 1·5 mm), the time taken to achieve
parameters in gas–solid reaction simulation have been 100% conversion has increased from 1·7 x 103 s to 4·5
considered. The isothermal model and kinetic x 103 s, as the radius of the grain sphere is increased
parameters similar to reduction of cobalt oxide (CoO) from 2·5 µm to 10 µm. This behaviour has been seen
with methane have been used for verifying parameters. in a number of reduction reactions. For example, the
The reduction of cobalt oxide with methane is as rate of reduction of nickel oxide by hydrogen is shown
follows: to increase due to grinding of nickel oxide.19 The shape
of grains also has been changed and its effect on the
CH4 + CoO → CO + 2 H2 + Co (35)
conversion–time curve was evaluated by the model.
We studied the kinetics of this reaction by Figure 2b shows the effect of grain shape factor
thermogravimetry, using synthetic porous pellets of variation for a slab-like pellet (Fp = 1). Spherical
pure CoO and various CH4 + Ar mixtures. The grains (Fg = 3) show higher conversion and shorter

a b

3 a) Effect of shape factor change, Fp = 1 (dashed line), Fp = 2 (dashed line), Fp = 3 (solid line) and (b) pellet
porosity change on the simulation shown for a slab-like pellet with spherical grains

Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 C15
Khoshandam et al. Simulation of non-catalytic gas–solid reactions

a b

c d

e f

g h i

4 Effect of changing initial porosity of pellet on the simulation (a), on gas reactant profile in the pellet, porosity in
the pellet, unreacted core radius and grain radius, for porosity = 40% (b, d, f and h), for porosity = 20% (c, e, g
and i)

C16 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114
 Khoshandam et al. Simulation of non-catalytic gas–solid reactions

time for completion of reaction than cylindrical (Fg = 2)

and slab (Fg = 1) grains and scales directly with surface
area/volume ratio.

Effect of pellet specifications

The effect of pellet shape factor on the conversion–time
curve has been evaluated by the model. Figure 3a shows
the results for spherical grains. As can be seen, a
spherical pellet made up of spherical grains (Fp = Fg = 3)
shows the highest conversion and shortest time for
completion of reaction but the difference with grains
taken as cylindrical (Fg = 2) and slab (Fg = 1) are
relatively small.
Indeed, porosity is one of the most important
specifications of the pellet. The effect of variation of
porosity on the conversion–time curve was also 5 Effect of changing Z in the simulation
investigated by the model. Figure 3b shows the effect
of changing porosity for a slab-like pellet (Fp = 1)
made up of spherical grains (Fg = 3). As the porosity grain size in the course of the reaction. Change in
of the pellet increased, the conversion of reaction has porosity of solid is indeed accompanied by a change in
been increased. Increasing the porosity also reduce the the effective diffusivity (see Equation 23).
completion time of the reaction. A 4-fold increase in The effect of changing the structure on
porosity from 20% to 80% cut down the 100% conversion–time behaviour has been studied using the
conversion time by a factor of 1.5-fold from 3·7 x 103 s model. Figure 5 shows the effect of different values of Z
to 2·4 x 103 s. Similar effects of porosity change on on the behaviour of conversion–time in a spherical pellet
conversion–time curves have been reported in the made up of spherical grains. With increasing Z values,
literature for reduction of nickel oxide by hydrogen.6 the time for 100% conversion increases. For very large Z
In the present model, as mentioned, effective values (> 2·2), conversion levels off with time before
diffusivity, conductivity and specific heat capacity of complete conversion. Such incomplete conversions have
the pellet are all functions of the porosity (Equations been experimentally observed for reaction of SO2 with
23, 28 and 27). Consequently, it is expected the limestone,20,21 and in the hydrofluorination of UO2.9 It
variation in porosity can strongly change the has been suggested that incomplete conversion is related
conversion–time curve. Figures 4a–4i illustrate the to ‘pore closure phenomenon’ whereby the porosity at
effect of initial porosity change on transient behaviour the surface of pellet becomes zero. As also seen in Figure
of gas concentration, local porosity, unreacted core 5, the fine grid should be used in the pellet in finite
radius and grain radius inside the pellet for a spherical volume method when Z is increased. The reason is to
pellet made up of spherical grains (Fp = Fg = 3). The achieve good solution accuracy in values of Z > 1.
concentration of reducing gas agent through the pellet
with high initial porosity (40%) reaches its maximum
value (equal to its value in bulk flow) quickly (Fig. VALIDATION OF THE MODEL
4b).On the other hand, the concentration in the pellet Comparison with analytical solutions – isothermal system
with low porosity (20%) changes slowly with time In this section, the mathematical model developed in
(Fig. 4c). Similar behaviour is also seen in the local this work is compared with analytical solutions given
porosity change through the pellet (Fig. 4d,e). by Wen1 for: (i) two limiting cases in gas–solid
Unreacted core radius of the grains in pellet surface reactions represented by the unreacted-core shrinking
layers shows a linear reduction due to time but in model and homogeneous model; and (ii) for the
deeper layers of the pellet, the variation is not linear intermediate condition between the two limiting
(Fig. 4f,g). Similar behaviour for changing local grain conditions. For a spherical pellet, assuming constant
radius are seen through the pellet (Fig. 4h,i). values for the diffusion coefficients, De and DAp, for a
first order reaction with respect to gas reactant A, the
Effect of structural change during reaction analytical solution under pseudo-steady state
If the solid product D occupies a larger volume than the approximation has been given by Wen.1 The solution
reactant B, a progressive decrease in porosity with was done in two stages: first, the period of reaction
conversion leads to cessation of reaction at incomplete prior to the formation of the product layer; and
conversion. The changes in grain dimensions determined second, the period following the product layer
by the parameter Z, defined as the ratio of the molar formation. If the resistance to the film diffusion
volume of the product to the reactant formulated as: around the pellet is negligibly small represented by a
dt B M D large Sherwood number, (Sh → ∞), for the first stage:
Z= (36)
bt D _1 - f D i M B X = 32 _ { coth { - 1i .i (37)
{
If Z < 1, the grains shrink during reaction;, if Z > 1,
swelling occurs and for Z = 1, there is no change in { = r p k vDC B0 (38)
e

Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 C17
Khoshandam et al. Simulation of non-catalytic gas–solid reactions

6 Comparison of conversion versus time; analytical 7 Comparison of conversion versus time; analytical
solution by Wen1 (unreacted-core shrinking model) solution by Wen1 (homogeneous model) and present
and present work work

i = k v C A0 t (39) prediction and analytical solution is shown in Figure 8.

In this case, agreement is adequate between model
where φ is Thiele modulus, θ is dimensionless time, kv
prediction and analytical solutions. Consequently, with
is reaction rate constant based on volume and CB0 is
attention to actual gas–solid reactions that cannot be
initial concentration of solid reactant.
completely described either by the unreacted-core
For the second stage:
shrinking model or by the homogeneous model,
3p c
X = 1 - p 3c + 2 9 {p c coth _ {p c i - 1C (40) prediction of the conversion by the present model is
{
expected to show good results for practical systems.
rcp
p c = rp (41)
Comparison with analytical solutions – non-isothermal
J N
i = 1 + c1 - D D e ln K p c sinh { O + system
Ap
m K sinh _ {p c i O For the analysis of non-isothermal systems, usually
L P there is no analytical solution, but under certain
{2 2
1 - p c i _1 + 2p c i + circumstances simple approximate solutions are
6 _
possible. For the case when the movement of the
De
D Ap _1 - p c i9 {p c coth _ {p c i - 1C (42) external surface toward the centre of the pellet is very
slow (concentration profile corresponds to the steady-
where ξc is dimensionless radius and rcp is radial state profile), under pseudo-steady state conditions
position of the reaction front in the pellet. the rate of diffusion of reactants across a boundary
The unreacted-core shrinking model can be used normal to the diffusion path is equal to the rate of
for porous pellets when the diffusion coefficient in the reaction within the boundary.
product layer is greater than the effective diffusion Consequently, by assuming pseudo-steady state
coefficient of gas through the pellet. Figure 6 shows conditions for the heat conduction in the pellet (in this
comparison between the analytical solution and
unreacted-core shrinking model prediction for the
case when DAp >> De (DAp/De = 100). There are
deviations between model prediction and exact
solution in conversions > 0·6.
The homogeneous model can also be used for
porous pellets when the pellet contains enough
voidage so that the gas reactant can diffuse freely into
the interior of the pellet. In this model, we assume the
effective diffusivity of the gas through the pellet is
equal to diffusion coefficient of the gas through the
product layer. Figure 7 compares the analytical
solution and model prediction for the case when DAp =
De. In this figure, deviations are observed between
conversion predicted by the model and analytical
solution for small Thiele modulus. Agreement is
acceptable for large values of Thiele modulus.
For conditions between two extreme models, the 8 Comparison of conversion versus time; analytical
comparison of conversion behaviour between model solution by Wen1 and present work

C18 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114
 Khoshandam et al. Simulation of non-catalytic gas–solid reactions

Table 8 Parameters used for the simulation of gas–solid

reaction in a non-isothermal system

Pellet Fp = 1, rp = 2·3 mm, ε0 = 70%,

De = 2·688 x 10–5 m2 s–1,
Ke = 0·1847 W m–1 K–1

Grain Fg = 3, r0 = 5 mm

Temperature Tb = 1073 K

Flow CAb = 0·0074 kmol m–3

Chemical reaction fG = k CA,

k = 1755·5 exp(–155560/RT),
∆H = 208 kJ mol–1

External heat and mass transfer Nu → inf

Sh → inf

9 Comparison of temperature profiles in a slab-like

pellet undergoing an endothermic reaction; analytical Table 9 Input parameters for carbon gasification reaction
solution (the mesh), and present work (the dots) with CO2

Fp 1
Fg 3
case the heat generated or consumed by the reaction is T (K) 1073–1373
all transferred across the boundary), the following P (atm) 1
analytical solution is obtained for the temperature r0 (mm) 3·75 x 10–5
rp (mm) 1·283–1·448
profile:5 ε0 66·5%
R V
S _ - DH i D e W Ke (W/m.K) 1·89
T - Ts = S Ke WW_ C As - C A i (43) CPe (kJ/kg.K) 1·49
S fu 0·8
T X ρB (kg/m3) 1800
Figure 9 compares the temperature profile between
MB (kg/kmol) 12
analytical solution and model prediction for a slab-
like pellet. The parameters are given in Table 8. The
temperature profile predicted by the model shows
good agreement with the profile calculated by source term in discretised equations as follows:
analytical solution.
# S U dV = SU DV = S u + S P C CO 2 (49)
CV

where
APPLICATION TO EXPERIMENTAL WORK df G
S u = f G (C -e C CO (50)
CO2 , old)
dC CO o
2
2, old

(C CO2, old)
Carbon gasification reaction with carbon dioxide
df G
The model used for gasification reaction with carbon Sp= e (51)
dC CO o 2
dioxide is: (C CO2, old)

When f G (C , old) is rate value due to _ C CO , old i,

CO2 2

C (solid) + CO 2 (gas) 2CO (gas) (44)

The experimental data reported by Abba and
Hastaoglu4 are given in Table 9. The reaction rate can
be described as:

fG = kPCO 2
(45)
1 + k1 PCO + k 2 PCO 2

where PCO and PCO are partial pressures of carbon

2

dioxide and carbon monoxide, respectively, and the

coefficients k, k1 and k2 are defined as:
k 8 mol kPa $ m 2 $ s B = 10 $ 8 exp _ - 219000 RT i
(46)

k18 kPa - 1 B = 3 $ 26 x 10 - 11 exp _ 253000 RT i (47)

k 2 8 kPa - 1 B = 1$ 74 x 10 - 3 exp _ 277000 RT i (48)

10 A comparison between simulation results with present
In this case, the source term in the transport program and experimental data for carbon gasification
Equations (2) and (12) is not first order with respect to reaction with CO2, the solid lines show simulation
gaseous reactant so we need to use a linear form of the results and circles show experimental data

Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 C19
Khoshandam et al. Simulation of non-catalytic gas–solid reactions

df G
e dC o is the
CO 2
(C CO2 , old)

value of the rate equation

derivative due to C CO , old and C CO , old is the concent-
2 2

ration of CO2 in the previous time step.

The comparison between simulation results and
experimental data is shown in Figure 10. There is
good agreement between mathematical modelling
results by the present program and experimental data
at different reaction temperatures.

Reduction of cobalt oxide with methane

The model has been applied to the reduction of cobalt
oxide with methane. Reduction of cobalt oxide with
methane can be described by the equation:
11 Application of the model for reduction of cobalt
CoO (solid) + CH 4 (gas) Co(solid) + CO (gas) + 2H 2 (gas)
oxide with methane
(52)
This reaction is endothermic (∆H = +208·4 kJ/mol)
and assumed to be first order with respect to methane complexity to numerical solutions as well. The
concentration. The kinetic data for this reaction were multiple reactions (between two solids and a gas) need
obtained using thermogravimetry techniques as more balance equations for solid reactants. This
follows:18 requires changes to parts of the program to get valid
results.
k = 1755 $ 5 exp _ - 155560/RT i (53) In situations where the gas–-solid reactions have
where k′ = 1755·5 m/s, Ea = 155 560 J/mol and R = complex mechanisms (for example in the cases where
8·314 J/mol K. Figure 11 compares the experimental solid reactant plays a catalytic role during the course
data and mathematical model solved with the finite of the reaction), detailed analysis of the reaction
volume technique in the present work for different mechanism and major modifications to the program
parameter v t , the generalised gas–solid reaction will be required. It should be mentioned that all non-
modulus. For the Equation (52), this parameter is catalytic gas–solid reactions can not be discussed with
defined as: the grain model and sometimes the reaction
mechanism is more complex than what has been
_1 - f 0 i kF p A g
t = Vp
v d Fg V g n (54) considered here.
Ap 2D e
Parameter v t incorporates both kinetic and structural
properties and is a measure of the relative magnitude CONCLUSIONS
of chemical reaction and diffusion rates. This In the present work, a mathematical model based on
parameter plays an important role in defining the the grain model is presented. The model uses transient
behaviour of gas–solid reactions; when v t (nominally state for mass and heat balance equations and also
t
v < 0·3), the overall rate of reduction is controlled by considers structural change such as porosity and grain
chemical kinetics (region 1), and when v t (nominally vt size changes in simulation of gas–solid reactions. The
> 2·0), the overall rate is controlled by product layer finite volume method was used for solving the
diffusion (region 3). Between these two extremes, both equations and a program in MATLAB was developed.
chemical kinetics and product layer diffusion are The effect of grain and pellet specifications, solid
significant (region 2).22 For reduction of cobalt oxide product porosity and reaction kinetics on the
with methane, there is good agreement between model gas–solid reactions were investigated using the
predictions and experimental data for conversions program.
below 0·9, for the values of v t in the range 0·2–0·6 Comparison between program results and
selected for the experimental conditions at 800–950°C. analytical solution, presented by Wen,1 shows that the
Often, the reaction can be more complex than model can be used for conditions between two limiting
considered here. In non-catalytic gas–solid reactions, conditions defined by unreacted core shrinking model
many examples can be found which contain reaction and homogeneous model. The program results also
kinetics of fractional order with respect to the gaseous compare favourably with analytical solution for
reactant (e.g. reduction of chromium oxide with pseudo steady state approximation for heat
hydrogen,23 or oxidation of zinc sulphide24). In these conduction in the pellet. Also, there is good agreement
systems, non-linearity can appear in the source term between simulation results and experimental data for
of transport equations. For such non-linear systems, the carbon gasification reaction with CO2 reported in
only numerical solutions may be used with correct the literature. The mathematical model was validated
linear assumptions for the source term in discretised with experimental data obtained in this work for
equations. Highly exothermic reactions can add reduction of cobalt oxide with methane and the results

C20 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114
 Khoshandam et al. Simulation of non-catalytic gas–solid reactions

show good agreement in regions where the overall rate introduction to computational fluid dynamics, the finite
of reaction is controlled by chemical kinetics or where volume method’, chapter 8, London, Adison Wesley
both chemical kinetics and product layer diffusion are Longman, 1995.
significant. 18. B. KHOSHANDAM: ‘Reduction of cobalt oxide, chromium
oxide, and manganese oxide with methane’, chapter 3,
Unpublished PhD Thesis, Amir-Kabir University of
Technology (Tehran Polytechnic), Tehran, Iran, 2003.
REFERENCES 19. H. Y. SOHN and M. E. WADSWORTH: ‘Rate processes of
1. C. Y. WEN: ‘Noncatalytic heterogeneous solid fluid extractive metallurgy’, New York, Plenum, 1979.
reaction models’, Ind. Eng. Chem., 1968, 60, 34–54. 20. M. HARTMAN and R. W. COUGHLIN: ‘Reaction of sulfur
2. P. A. RAMACHANDRAN and L. K. DORAISWAMY: dioxide with limestone and the influence of pore structure’,
‘Modeling of noncatalytic gas-solid reactions’, AIChE J., Ind. Eng. Chem. Proc. Des. Dev., 1974, 13, 248–253.
1982, 28, 881–900. 21. M. HARTMAN and R. W. COUGHLIN: ‘Reaction of sulfur
3. G. UHDE and U. HOFFMANN: ‘Noncatalytic gas–solid dioxide with limestone and the grain model’, AIChE J.,
reactions: modelling of simultaneous reaction and 1976, 22, 490–498.
formation of surface with a nonisothermal crackling core 22. J. SZEKELY, C. I. LIN and H. Y. SOHN: ‘A structural model
model’, Chem. Eng. Sci., 1997, 52, 1045–1054. for gas-solid reactions with a moving boundary. V An
4. I. A. ABBA and M. A. HASTAOGLU: ‘Modelling of multi experimental study of the reduction of porous nickel-oxide
gas-solid reactions: effect of inert solids’, Trans IChemE, pellets with hydrogen’, Chem. Eng. Sci., 1973, 28,
1997, 75 Part A, 33–41. 1975–1989.
5. J. SZEKELY, J. W. EVANS and H. Y. SOHN: ‘Gas–solid 23. W. F. CHU and A. RAHMEL: ‘The kinetics of reduction of
reactions’, Chapter 4, New York, Academic Press, 1976. chromium oxide by hydrogen’, Met. Trans. B., 1979, 10B,
6. J. SZEKELY and J. W. EVANS: ‘Studies in gas-solid 401407.
reactions: part II. An experimental study of nickel oxide 24. K. J. CANNON and K. J. DENBIGH: ‘Studies on gas-solid
reduction with hydrogen’, Met. Trans., 1971, 2, 1699–1710. reactions. I The oxidation of zinc sulfide’, Chem. Eng. Sci.,
7. J. SZEKELY and C. I. LIN: ‘The reduction of nickel oxide 1957, 6, 145–154.
discs with carbon monoxide’, Met. Trans., 1976, 7B, 493.
8. J. W. EVANS, S. SOHN and C. E. LEON-SUCRE: ‘The
kinetics of nickel oxide reduction by hydrogen,
measurements in a fluidized bed and in a gravimetric
Authors
apparatus’, Met. Trans., 1976, 7B, 55–65. B. Khoshandam has a PhD in chemical engineering in 2004
9. E. C. COSTA and J. M. SMITH: ‘Kinetics of noncatalytic, from Amir-Kabir University (Tehran Polytechnic). He has a
nonisothermal, gas-solid reactions: hydrofluorination of BS degree in chemical engineering – oil industry process
uranium dioxide’, AIChE J., 1971, 17, 947–958. design in 1992 from Iran University of Science and
10. J. B. GIBSON and D. P. HARRISON: ‘The reaction between Technology and MS dgree in chemical engineering in 1997
hydrogen sulfide and spherical pellets of zinc oxide’, I & from Amir-Kabir University. He is a researcher in gas-solid
EC, Proc. Des. Dev., 1980, 19, 231. reactions and a programmer in MATLAB.
11. F. PATISSON, M. G. FRANCOIS and D. ABLITZER: ‘A non-
isothermal, non-equimolar transient kinetic model for gas- R. V. Kumar is a senior lecturer in the Department of
solid reactions’, Chem. Eng. Sci., 1998, 53, 697–708. Materials Science and Metallurgy and a Fellow of Trinity
12. F. LIU and S. K. BHATIA: ‘Solution techniques for Hall at University of Cambridge, UK. He obtained hid PhD
transport problems involving steep concentration from McMaster University in Canada. His research
gradients: application to noncatalytic fluid solid reactions’, interests include thermodynamics and kinetics of
Comput. Chem. Eng., 2001, 25, 1159–1168. metallurgical and materials processing, metal refining and
13. N. WAKAO and J. M. SMITH: ‘Diffusion in catalyst pellets’, recycling, electrochemistry, solid state ionic sensors for
Chem. Eng. Sci., 1962, 17, 825–834. process and pollution control, processing of sensors and
14. E. A. MASON, A. P. MALINAUSKAS and R. B. EVANS: actuators, and sustainable technologies.
‘Flow and diffusion of gases in porous media’, J. Chem.
Esmail Jamshidi is Professor of the Department of Chemical
Phys., 1967, 46, 3199–3216.
Engineering, Amir-Kabir University (Tehran Polytecgnic).
15. R. B. BIRD, W. E. STEWART and E. N. LIGHTFOOT:
He has a BS degree in chemical and petroleum engineering
‘Transport phenomena’, chapter 16, New York, Wiley,
fron Tehran University in 1964, MS and PhD degree in
1960.
chemical engineering and materials science from Wayne
16. W. E. RANZ and W. R. MARSHALL: ‘Evaporation from
State University in 1971 and 1973 respectively. His research
drops’, Chem. Eng. Prog., 1952, 48, 141–146, 173–180.
interests are methane reactions wit inorganic chemicals.
17. H. K. VERSTEEG and W. MALALASEKERA: ‘An

Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114 C21
Khoshandam et al. Simulation of non-catalytic gas–solid reactions

APPENDIX – NOTATION

A, C Gaseous reactant and product, respectively r Radial position in the grain (m)
Ag, Ap External surface area of individual grain and rc Radial position of the reaction front in grain (m)
pellet, respectively (m2) rcp Radial position of the reaction front in
B, D Solid reactant and product, respectively analytical solution (m)
b, c, d Stoichiometric coefficients rg Radius of the grain (m)
CB0 Initial concentration of solid reactant rp Pellet radius for spherical and cylindrical pellets
(kmol m–3) and half-thickness for slab-like pellets (m)
CA Gaseous reactant concentration in the pellet S Average value of source S over the control
(kmol m–3) volume
CAb, CAs Gaseous reactant concentration in the bulk of SΦ Source term in Equation (31)
gas and on the external surface of pellet, Su, Sp Constant coefficients in the linear source model,
respectively (kmol m–3) Equation (32)
C CO2 Concentration of carbon dioxide in the bulk of T Temperature (K)
the gas phase t Time (s)
C CO , old Concentration of carbon dioxide in the bulk of
2 u∞ Velocity of gas stream (m s–1)
the gas phase in the previous time step Vg, Vp Volume of individual grain and pellet,
CPe Effective specific heat capacity of the pellet respectively (m3)
(J kmol–1 K–1) X Overall solid conversion for the pellet
CPf Specific heat capacity of gas stream x Local conversion in the pellet
(J kmol–1 K–1) Z Volume of solid product from unit volume of
CPs, CPg Specific heat capacity of solid and gas, solid reactant
respectively (J kmol–1 K–1) z Position in the pellet, radial position in spherical
CV Control volume and cylindrical pellets and axial position in
De Effective diffusion coefficient of gaseous slab-like pellets
reactant through the pellet (m2 s–1)
DAp Diffusivity of gaseous reactant A in the product
GREEK LETTERS
layer (m2 s–1)
Γ Diffusion coefficient of property Φ in Equation (31)
DAK, DAM Knudsen and molecular diffusivity of gaseous
∆H Heat of reaction (J kmol–1)
reactant A through the pellet, respectively
ε Porosity of the pellet
(m2 s–1)
fu Emissivity of the pellet
dp Diameter of pellet (m)
θ Dimensionless time
Ea Activation energy (J kmol–1)
µf Viscosity of gas stream (kg m–1 s–1)
Fp, Fg Pellet and grain shape factor, respectively
ξc Dimensionless radius defined by Equation (41)
fG Rate of surface reaction (kmol m–2 s–1)
ρ Density (kg m–3)
hc Convective heat transfer coefficient (W m–2 K–1)
ρf Density of gas stream (kg m–3)
Ke Effective thermal conductivity of the pellet
σ Stefan-Boltzmann constant 5·67 10–8 (W m–2 K–4)
(W m–1 K–1) t
v Generalised gas–solid reaction modulus defined
Kf Thermal conductivity of gas stream (W m–1 K–1)
by Equation (46)
Ks, Kg Thermal conductivity of solid and gas,
τ Tortuosity coefficient of the pellet
respectively (W m–1 K–1)
Φ A scalar property in Equation (31)
k Reaction rate constant (m3n-3m-2 kmol1-m-n s–1)
ΦP Property Φ at a general nodal point P
kv Reaction rate constant based on volume
ϕ Thiele modulus
(m3(m+n-1) kmol1-m-n s–1)
kc Convective mass transfer coefficient (m s–1)
k0 A parameter defined in Equation (25) DIMENSIONLESS NUMBERS
k1 A coefficient in Equations (21) and (45) Nu Nusselt number = hc dp/Kf
k2 A coefficient in Equation (45) Pr Prandtl number = CPf µf/Kf
k′ Frequency factor (kmol1-m-n s–1 m3n-3m-2) Re Reynolds number = ρf u∞ dp/µf
M Molecular weight (kg kmol–1) Sc Schmidt number = µf/ρf DAM
n, m Exponent of the reaction rate dependence on Sh Sherwood number = kc dp/DAM
the reactant concentrations
nd Number of solid grains per unit volume of SUBSCRIPTS
porous solid (m–3) 0 Initial conditions
P Total pressure b Bulk of gas flow
PCO , PCO Partial pressures of carbon dioxide and
2 i Reaction front condition
carbon monoxide s External surface of the pellet
R Ideal gas constant (J mol–1 K–1) w Inner walls of the reactor

C22 Mineral Processing and Extractive Metallurgy (Trans. Inst. Min. Metall. C) March 2005 Vol. 114