New modelling approach to liquid–solid reaction kinetics: From ideal particles to real particles

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

Contents lists available at ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

New modelling approach to liquid–solid reaction kinetics: From ideal particles to real particles Tapio Salmi ∗ , Henrik Grénman, Johan Wärnå, Dmitry Yu. Murzin Laboratory of Industrial Chemistry and Reaction Engineering, Process Chemistry Centre, Åbo Akademi, Biskopsgatan 8, FI-20500 Turku-Abo, Finland

a b s t r a c t Typical features of liquid–solid reactions were reviewed: reaction kinetics, mass transfer effects and particle morphology. It was concluded that classical liquid–solid models based on ideal, non-porous geometries (sphere, infinite cylinder, slab) cannot satisfactorily describe real reactive solid particles with various surface defects, such as cracks, craters and limited porosity. Typically a too low reaction order for the reactive solid is predicted by the classical models. The surface morphology can be revealed by electron microscopy, which gives inspiration to develop new mathematical models for reactive solids. Generalized product layer and shrinking particle models were developed for arbitrary geometries deviating dramatically from ideal geometries. The approach is based on the use of a shape factor (a), which can obtain high values for particles with surface defects. The modelling approach was illustrated with numerical simulations and an experimental case, leaching of zinc sulphide with ferric ions. A general methodology for investigation of fluid–solid reactions was proposed. © 2013 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. Keywords: Fluid–solid reactions; Kinetics; Diffusion; Modelling; Morphology; Shape factor; Particle size distribution

1.

Introduction

1.1.

Solid–fluid reactions in general

The interest for gas–solid reactions with reactive solids has been large in the recent years. One of the most important applications is the combustion of solid fuels; thus a lot of mathematical modelling work has concerned the interaction between gases and reactive solids. Good reviews on models of gas–solid reactions are given, for instance, by Szekely et al. (1976) and Trambouze et al. (1988). Classical models for gas–liquid systems are shrinking particle and product layer models, for which ideal particle shapes (non-porous spheres, long cylinders and flakes) are assumed. The research effort devoted to modelling of solid–liquid systems has been much more modest, even though numerous solid–liquid processes are applied in practise. Just to mention few examples, solid–liquid reactions appear in gastronomy, in metallurgy, in the production of organic fine chemicals as well as in the

∗

synthesis of pharmaceuticals. Many metals are extracted from ores by hydrometallurgical processes, in which the metals are selectively dissolved (Habashi, 2005). Chemical Kraft pulping (Christensen et al., 1983; Santos et al., 1997; Andersson, 2003; Yang and Liu, 2005; Borgen et al., 2007, 2009; Grénman et al., 2010) is used to produce paper by selectively extracting lignin from wood chips; cellulose fibres remain in solid phase and are processed further for papermaking (RISI Report, 2009). Fragrances, pharmaceuticals and food ingredients are extracted from plants. Very recently, hot-water extraction of valuable components, such as hemicelluloses from wood chips has got increased attention (Willför et al., 2004, 2005; Grénman et al., 2011). The solid–liquid interaction – intrinsic kinetics coupled to mass and heat transfer effects – plays a central role in all these processes. In general, it can be stated that similar theories can be applied to both gas–solid and liquid–solid processes. Therefore, they are called fluid–solid reactions in the sequel.

Corresponding author. Tel.: +358 2 2154427; fax: +358 2 2154479. E-mail addresses: Tapio.Salmi@abo.fi, tsalmi@abo.fi (T. Salmi). Received 13 April 2013; Received in revised form 30 July 2013; Accepted 2 August 2013 0263-8762/$ – see front matter © 2013 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cherd.2013.08.004

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

Experimental activities and mathematical modelling of fluid–solid reactions follow the general principles of chemical engineering. The starting point in designing a process involving chemical transformations should be based on a detailed understanding of the reaction thermodynamics and kinetics. Thermodynamics limits the extent of any chemical reaction by setting the boundaries of chemical equilibria. The major factors governing the rates of chemical reactions are the concentrations or activities of the reacting species, temperature, pressure as well as mass and heat transfer limitations. The investigation of reaction rates is usually performed in laboratory reactors due to economic reasons, but particularly because well-controlled reaction conditions can be obtained in a small scale (constant temperature and pressure, gradientfree conditions, rapid mass and heat transfer). In industrial scale, however, mass and heat transfer rates influence the overall reaction rates. Mass transfer effects are more dominant in solid–liquid reactions compared to solid–gas reactions as the mass transfer is slower and the heat transfer is faster in liquid phase than in the gas phase. A common methodology in investigating the kinetics of chemical reactions is to carry out experiments in a batch, semi-batch or continuous reactors. The progress of the reaction is often measured by observing the change in the weight of the solid phase or by measuring the product concentrations in the liquid phase. Product identification from the fluid phase samples can be performed with various analytical techniques, such as gas and liquid chromatography, titrimetry, flow injection analysis, atomic absorption spectrometry (AAS) and inductively coupled plasma (ICP). The influence of external mass transfer can be experimentally determined by, for example, performing experiments at different stirring rates until no influence on the reaction rate can be obtained by further increasing the agitation. Internal mass transfer is relevant only for porous solids and this can be experimentally investigated, for example, by performing experiments with different particle sizes. A reactor system for precise studies of liquid–solid kinetics has recently been proposed by Grénman et al. (2011). The system consists of stirred tanks in series and the liquid flow is recycled with high flow rates. Each tank is perfectly backmixed and has the same liquid phase composition, because the recycle flow rate is very high compared to the reaction rate. At the first time moment, the last reactor is decoupled from the system and both solid and liquid phase are analyzed. At the second sampling time moment the following reactor is decoupled etc. The reactor system is displayed in Fig. 1.

1.2.

Overview of fluid–solid reaction kinetics

1.2.1.

Specific feature: Solid surface area

There is a fundamental difference between homogeneous and fluid–solid kinetics: The reaction rate of a solid-phase compound dispersed in the fluid does not directly depend on the concentration of the solid in the fluid phase, but rather on the reactive surface area. This fact complicates the work very much, since the quantification of the surface area is not as straightforward as the concentration measurements. The surface area does not usually decrease linearly with the reactant conversion as the concentrations do, which makes the kinetics to behave differently and complicates the interpretation of experimental data. Furthermore, the reactive, accessible surface area is not always equal to the total surface area, e.g. certain crystal sites might react differently. The solid

1877

Fig. 1 – Stainless steel cascade reactor system constructed for enhanced wood extraction studies. The five 200 mL Parr reactors were heated with heating jackets and heating cables. The liquid was circulated with a powerful membrane pump. component may also be a heterogeneous one comprising of a number of compounds as is often the case in, for example, hydrometallurgy and pulping. Even phase changes in the solid may sometimes occur due to extreme reaction conditions. These factors are, however, usually proportional to the total surface area. The specific surface area of solids (, in m2 /g solid) in a static state can be measured with nitrogen physisorption and the adsorption isotherms can be interpreted e.g. with BET and Dubinin theories. This approach is rather time consuming which in practice restricts its use; typically the surface area is measured before and after the kinetic experiments. The accuracy of surface area measurements is not the best possible for low areas typical for fluid–solid systems. The morphology of the solid phase can be altered in the sample preparation as the samples usually have to be dried prior to analysis. The particle size distribution of a solid can be studied in situ with laser diffraction or dynamic light scattering. Theoretical functions, such as gamma function can sometimes be fitted to experimentally recorded particle size distributions. A typical example is shown in Fig. 2. If the morphology of the particles is known to a reasonable extent, the surface area can be estimated by varying the particle size. No more accurate method exists for this purpose and values obtained can be very imprecise especially for particles for which the porosity changes during the reaction. Overall, direct in-situ measurements of surface properties of micrometer scale particles during a reaction remain a great challenge.

1.2.2.

Modelling aspects—SEM reveals the truth

Certain assumptions need to be made about the behaviour of the solid for modelling purposes, if the surface area during the reaction cannot be quantified directly, which is the typical case. The most commonly used assumptions are based on relating the change in surface area to the reactant conversion, which can be more easily measured than the surface area. This correlation is based on the morphological changes of the solid during the reaction. The particle morphology is typically assumed to be uniform and it is assumed to be of some ideal, non-porous shape, e.g. sphere, cylinder or slab. Based on elementary geometry, it can be shown that the surface area of

1878

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

10 9

Frequency (counts/s)

8 7 6 5 4 3 2 1 0 0

50

100

150

200

250

Diameter ( m)

Fig. 2 – The fit of the gamma function to the experimental particle size distribution of gibbsite. spheres (4R2 ) changes with conversion, i.e. the volume of the spheres (4R3 /3), in a certain way, independent of the particle size. For the solid structure, shrinking particle, product layer (shrinking core) and grain models have been proposed (see e.g. Trambouze et al., 1988 and Levenspiel, 1999). According to the shrinking particle model, all products dissolve to the fluid phase or crack away from the surface as soon it is formed. Product layer model assumes that the solid product forms a porous layer around the solid particle. This layer causes increasing mass transfer resistance as the reaction progresses, since the reactant molecule has to diffuse through this layer before reaching the reaction surface. The grain model assumes a raspberry structure of the solid phase: The fluid molecules diffuse through the structure to non-porous microparticles, where the reactions take place. The shrinking core (product layer) and shrinking particle models are illustrated in Fig. 3 (Debasish and Bandyopadhyay, 2011; Hossein et al., 2010; Suchun and Nicol, 2010; Fabiano et al., 2010). If the whole particle is porous from the beginning, a porous particle model appears, in analogy with the treatment of the reaction–diffusion problem in heterogeneous catalysis. This model is rather seldom used for fluid–solid reactions, because it is generally agreed that a steep reaction front appears inside the solid particle, for instance, in combustion processes. However, when solids react with liquids at moderate temperatures, such as in the leaching of biomass or delignification of wood, the porous particle model becomes relevant; the change of the

porosity is the crucial issue, since solid material is dissolved out from the porous particle. For different particle geometries, mathematical relations between the reaction time, surface area and the solid reactant conversion (˛ = conversion of the solid substance) can be derived. Some examples are collected in Table 1. The expressions are valid for isothermal (semi)batch reactors; the fluid phase concentrations are assumed to be constant, which is a typical case when a thermogravimetric equipment is used: an excess of gas flows through a measurement chamber, where the solid sample is placed. The models listed in Table 1 can be tested by plotting experimental data (factor g) versus the reaction time. If the plot is a straight line, it is usually concluded that the model is valid. The main weakness in this methodology is that the solid particles seldom follow ideal morphology. Moreover, the discrimination of the models can be challenging as several models can fit the experimental data reasonably well, as the parameter values are adjusted, especially in case of considerable experimental scattering. A common feature for the models proposed so far is the hypothesis of ideal geometry, i.e. slab, infinitely long cylinder and spherical particle is assumed. For many cases, this is not the truth. A new approach is needed: If the experimental case does not fit to Table 1, what to do? The main reason for the deviation from the ideality (Table 1) is often the highly non-ideal structure of the solid material: modern microscopic techniques, particularly scanning electron microscopy (SEM) have revealed that real particles often considerably differ from the ideal ones. Surface defects, such as cracks and craters, as well as porosity affect the particle structure and the observed reaction kinetics. Thus, there is a risk that assuming an ideal geometry for a solid particle, a wrong reaction order with respect to the solid substance is predicted. For example, the hypothesis of a perfect, non-porous spherical particle free from any defects implies the reaction order 2/3 with respect to the solid component (Table 1, case 5), but a higher reaction order is expected, if the surface is rough and/or the particle contains some pores (the mathematical treatment is provided in Section 5). Some examples of real, non-ideal surfaces are displayed in Fig. 4. The SEM images strongly suggest that the particles deviate from the well-known ideal geometries. The measurements of specific surface areas by gas adsorption (physisorption) have also revealed that the surface areas of reactive solids can clearly exceed the areas predicted by ideal surfaces. For instance, a surface area of 2 m2 /g – a superficially low one – has been measured by nitrogen adsorption on solid, spherical particles. If the density of the solid material is around 1500 kg/m3 and a perfect, non-porous structure is presumed,

Fig. 3 – Shrinking core (product layer) and shrinking particle models. In shrinking particle model, no product layer appears.

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

1879

Table 1 – Factors for solid particles (Dickinsin and Heal, 1999; Órfão and Martins, 2002). g(˛)

f(cS )

1 2 3 4 5 6 7 8 9

−ln(1 − ˛) (1 − ˛)−1/2 − 1 (1 − ˛)−1 − 1 1 − (1 − ˛)1/2 1 − (1 − ˛)1/3 1 − (1 − ˛)2/3 [1 − (1 − ˛)1/3 ]2 1 − 2˛/3 − (1 − ˛)2/3 [1/(1 − ˛)1/3 − 1]2

cS /c0S (cS /c0S )3/2 (cS /c0S )2 (cS /c0S )1/2 (cS /c0S )2/3 (cS /c0S )1/3 (cS /c0S )2/3 /(1 − (cS /c0S )1/3 ) (cS /c0S )1/3 /(1 − (cS /c0S )1/3 ) (cS /c0S )5/3 /(1 − (cS /c0S )1/3 )

10 11 12

[1 − (1 − ˛)1/2 ]2 1/(1 − ˛)1/3 − 1 1 − 3(1 − ˛)2/3 + 2(1 − ˛)

(cS /c0S )1/2 /(1 − (cS /c0S )1/2 ) (cS /c0S )4/3 (cS /c0S )1/3 /(1 − (cS /c0S )1/3 )

Type of model First-order kinetics Three-halves-order kinetics Second-order kinetics One-half-order kinetics; two-dimensional advance of the reaction interface Two-thirds-order kinetics; three-dimensional advance of the reaction interface One-thirds-order kinetics; film diffusion Jander; three-dimensional Crank–Ginstling–Brounshtein, mass transfer across a nonporous product layer Zhuravlev–Lesokhin–Tempelman, diffusion, concentration of penetrating species varies with ˛ Jander; cylindrical diffusion Dickinson, Heat transfer across the contacting area Shrinking core, product layer (different form of Crank–Ginstling–Brounshtein)

Factor g(˛) (˛ = conversion of the solid = 1 − cS /c0S ) is directly proportional to the reaction time in an isothermal batch reactor, in which the liquid or gas phase component concentration is assumed to be constant. Factor f(cS ) explains the reaction order with respect to the solid substance (S) concentration. Example: dcS /dt = −kcS (case 1) gives after integration −ln(cS /c0S ) = kt, where cS /c0S = 1 − ˛

the average particle radius calculated from these data is 1 ␮m only! Typically the particle radii observed by microscopy are much higher, in the range of 10–1000 ␮m or even more. The specific surface area measured by gas adsorption and the particle size observed by microscopy can be explained only by surface defects and some degree of porosity. Here we present mathematical models for general particle geometries, which allow the description of ideal non-porous particles and completely porous particles, and all intermediates between these extreme cases. The impact of this generalized model for shrinking particle and product layer models is discussed thoroughly and the modelling approach is illustrated with some common rate equations applied to fluid–solid kinetics.

2.

Generalized product layer model

2.1.

General characteristics of the product layer model

A generalized product layer model is presented here. The model is applicable to any geometry, including a non-ideal one, rough surfaces with cracks, craters and pores. The generalization is based on the use of a shape factor (a), which gives a relation between the accessible surface area to the particle volume. For an ideal, non-porous sphere (radius = R), the

shape factor is a = (Ap /Vp )R, where Ap is the outer surface area (Ap = 4R2 ) and Vp is the particle volume (Vp = (4/3)R3 ). Thus the shape factor becomes a = 3. It can be intuitively understood, that the shape factor can obtain much higher values for non-ideal particles (see Fig. 4), as the outer surface area can be much higher than predicted from the simple geometry. Based on this principle, model equations will be derived for solid particles. For the sake of illustration, a uniform particle size is presumed, but the model is easily extended to particle size distributions.

2.2. Diffusion through porous product layer in an non-ideal particle For the sake of illustration, an isotropic product layer is assumed. The mass balance of a diffusing fluid component in the product layer is described by the mass balance equation, (Ni A)in = (Ni A)out +

dni dt

(1)

where Ni is the diffusion flux and A is the surface area available for diffusion transport. The infinitesimal difference between the fluxes (out–in) is denoted by (Ni A) and the amount of substance is ni = cLi V, where V is the volume of the volume

Fig. 4 – Real solid particle: SEM image of zinc sulphide (ZnS) after some leaching with ferric ions (Fe3+ ).

1880

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

element. The volume element is let to approach zero, after which the balance equation becomes dcLi d(Ni A) =− dt dV

(2)

The further treatment depends on the detailed geometry and the diffusion model applied. If the simplest model, the random pore model with the concept of an effective diffusion coefficient for each component (Dei ) is applied, the flux is given by Fick’s law, Ni = −Dei

dcLi dr

(3)

where r is the length coordinate. For an ideal spherical geometry, the area and volume elements are 4r2 and 4r2 dr, respectively. For an infinitely long cylinder, these elements are 2rL and 2rLdr, respectively. For a slab, the corresponding elements are L2 and L2 dr. The flux expression (3) is inserted in the balance equation (2). It can easily be shown that the balance equation can be written in a general form dc d(Dei (dcLi /dr)ra−1 ) = Li dt ra−1 dr

(4)

where the shape factor is a = 1 for a slab, a = 2 for a cylinder and a = 3 for a sphere. Furthermore, the shape factor represents the ratio between the outer surface area (AP ) and the volume of the particle (VP ): a=(

AP )R VP

(5)

where R is the characteristic dimension, the original radius of the unreacted particle; for instance the radius of the cylinder or the sphere. For ideal surfaces, the value of the shape factor is maximally a = 3 (ideal spheres), but for non-ideal particles, this value can be exceeded due to surface defects and porosity. Pseudo-steady state conditions are usually assumed for product layers and thin films around the product layers. If the diffusion coefficient is assumed to be locally constant, the balance equation (4) is simplified to

 Dei

d2 cLi dr2

+

(a − 1) dcLi r dr

 =0

(7)

(8)

where C is the second integration constant. The further treatment of Eq. (8) depends on the boundary conditions applied. At the outer surface of the particle, the concentration is denoted ∗ (r = R). Thus we have according to Eq. (7), by cLi ∗ dcLi

dr

= CR1−a

(9)

(10)

giving the boundary condition



b kLi cLi −

C R2−a − C 2−a



= Dei CR1−a

at r = R

(11)

On the other hand, the solution (8) should be valid at the s , reaction surface (r), where the concentration is denoted by cLi s cLi =

C 2−a + C r 2−a

(12)

The integration constant C is eliminated from Eqs. (11) and (12), and the integration constant C is solved: C=

b s cLi − cLi

(13)

(R2−a − r2−a )/(2 − a) + Dei R1−a /kLi

This expression is inserted in Eq. (7) to obtain an equation for the component flux: Ni =

b s −Dei (cLi − cLi )

r((r/R)

a−2

− 1)/(2 − a) + Dei (r/R)

a−1

/kLi

(14)

After introducing the Biot number (BiMi ) for mass transfer, BiMi = (kLi R)/Dei and rearranging Eq. (14), a general flux expression is obtained: Ni =

b −(a − 2)Dei (cLi − cisLi )

R(1 − (1 − (a − 2)/BiMi )(r/R)

a−2

)(r/R)

(15)

It should be noticed, that Eq. (15) is valid for an arbitrary geometry, even for non-integer values of the shape factor (a). The models of ideal particles are obtained as special cases of Eq. (15). For a spherical particle, a = 3, and Eq. (15) becomes

(6)

where C is the first integration constant. A second integration gives the general solution C 2−a cLi = + C r 2−a

b ∗ Ni = −Dei CR1−a = −kLi (cLi − cLi )

Ni =

This balance equation can be solved for an arbitrary geometry (non-integer values of the shape factor) by integrating it twice. After the first integration, the concentration gradient is obtained, dcLi = Cr1−a dr

On the other hand, the fluid film resistance prevails at the outer surface of the particle. The fluxes expressed in different ways are equal,

b s −Dei (cLi − cLi )

(16)

R(1 − (1 − 1/BiMi )(r/R))(r/R)

For an infinitely long cylinder, a = 2, and Eq. (15) cannot be used as such. We denote a − 2 = z and consider the function F(z) = f (z)/g(z) =

z 1 − (1 − z/BiMi )(r/R)

z

(17)

F(z) → 0 as z → 0. Thus the rule of l’Hopital is applied: the limit of F is f (z)/g (z). It can be shown that f (z) = 1 and g (z) = 1/BiMi − ln(r/R). Thus the flux expression for an infinitely long cylinder becomes Ni =

b s −Dei (cLi − cLi )

R(ln (r/R)

−1

+ 1/BiMi )(r/R)

(18)

The same solution can be obtained by inserting a = 2 into Eq. (6) and solving it analytically. For an ideal slab (a = 1), Eq. (15) gives after some rearrangement the simple solution Ni =

s −Dei (cLbLi − cLi )

R(1 + 1/BiMi − r/R)

(19)

1881

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

It should be noticed that rapid mass transfer through the fluid film surrounding the particle implies that kLi → ∞, BiMi → ∞ and 1/BiMi → 0 in Eqs. (15)–(19).

2.3.

=

Interaction of surface reaction and diffusion

At the reaction surface (at the coordinate position r), the following balance equation can be written for each fluid-phase component:



By using a dimensionless number, a generalized Thiele modulus, −R

−1 = s ik Rk (cLi )

b s (a − 2)Dei (cLi − cLi )

R(1 − (1 − (a − 2)/BiMi )(r/R)

a−2

)(r/R)

+

S 

s ik Rk (cLi )=0

(21)

k=1

s The concentrations (cLi ) at the surface can be solved from the above equation by inserting the rate expressions (Rk ). For simple cases of reaction kinetics, such as zero, first and second order kinetics, an analytical solution of Eq. (21) is easily obtainable. The dimensionless radius (r/R) can be related to the amount of substance of the solid material. If we denote the initial amount of the solid material by n0j and the amount of unreacted solid at an arbitrary time (t) by nj , these quantities correspond to the particle volumes accordingly,

=

Vj

(22)

V0j

For a spherical particle is valid Vj /V0j = (r/R)3 , for a long cylinder Vj /V0j = (r/R)2 . For a general geometry, the relation Vj /V0j = (r/R)a is valid, which gives finally r = R



nj

a−2 (1 − (1 − (a − 2)/BiMi )(nj /n0j )

1−2x

)(nj /n0j )

x

(27)

the following expression is obtained,

where Rk denotes the chemical reaction rates. The reaction surface is assumed to be so thin that no accumulation of substances takes place; thus a pseudo-steady state hypothesis is applicable. After inserting the general flux expression into Eq. (15), the following relation is obtained between the chemical rates and the diffusion flux:

nj

(26)

(20)

k=1

n0j

 R (1) k=1 ik k

b Dei cLi

where Rk (1) is the rate at yi = 1, and a parameter

S

Ni A =

S

1/a (23)

n0j

S Sk=1

1 − yi − 

 R (1) k=1 ik k

s − cLi ) 1−2/a 1/a R(1 − (1 − (a − 2)/BiMi )(nj /n0j ) )(nj /n0j )

+

S 

1 − yi −

 R(y) R(1)

3.

Generalized shrinking particle model

3.1. model

General characteristics of the shrinking particle

3.2.

By denoting 1/a = x and introducing a dimensionless cons b /cLi , Eq. (24) becomes centration yi = cLi (a − 2)(1 − yi ) (1 − (1 − (a − 2)/BiMi )(nj /n0j )

1−2x

)(nj /n0j )

x

+

R

S

 R (y ) k=1 ik k i b Dei cLi

=0

=0

(29)

Interaction of chemical reaction and diffusion

According to the principles described above, the balance equation for each fluid-phase component at the reaction surface can be written as follows, provided that the mass transfer through the fluid film can be described with Fick’s law: S 

s ik Rk (cLi )=0

(30)

k=1

s ik Rk (cLi )=0

(24)

(28)

For the shrinking particle model, the treatment is much simpler, since the product layer is missing—thus the chemical reactions and film transfer resistance remain only. A general geometry is assumed here, with all values of the shape factor (a) possible. The balance equation can be written analogously with the balance for the product layer model, i.e. what is transported to the reaction surface is consumed by the reactants. For the reaction products, the opposite reasoning is valid. A pseudo-steady state is assumed on the reaction surface.

b s kLi (cLi − cLi )+

k=1

=0

from which the unknowns, the dimensionless concentrations (yi ) are solved iteratively or analytically. For the case that only one reaction proceeds in the system, Eq. (28) is simplified to

The algebraic equation system to be solved with respect to each fluid-phase component becomes now b (a − 2)Dei (cLi

ik Rk (y)

The mathematical analogy with the product layer model is evident; Eq. (30) can be brought to a dimensionless form

S

 R (y) k=1 ik k b kLi cLi

1 − yi +

=0

(31)

s b where yi = cLi /cLi . Furthermore, we can write

S

(25) 1 − yi −



Sk=1

ik Rk (y)

 R (1) k=1 ik k

=0

(32)

1882

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

analogously to Eq. (28). The dimensionless Thiele modulus is defined by

S

=

−

 R (1) k=1 ik k b kLi cLi

(33)

For the case that only one chemical reaction takes place in the system, Eq. (32) is simplified to  R(y) R(1)

1 − yi −

=0

(34)

The equation system (32) is solved iteratively or analytis cally with respect to each concentration cLi . Some specific examples will be considered in next chapter.

4. Analytical solution of the models for some specific kinetics

s2 kcLA s 1 + ˇcLA

(35)

s2 kcLA

(36)

s kcLA R = s 1 + ˇcLA

(37)

s R = kcLA

(38)

R = k

(39)



R =

The method is principally very straightforward. For instance, the rate equation (35) is inserted in the general balance Eq. (29) s b giving (yi = cLi /cLi ) b )y2  (1 + ˇcLA b 1 + ˇcLA y

=0

(40)

b The dimensionless quantities ˇ = ˇcLA and ˛ =  (1 + are introduced and the second-degree equation is solved. The physically meaningful solution of Eq. (40) is

ˇ )

y=

2



2

(1 − ˇ ) + 4(˛ + ˇ ) + 1 − ˇ

(41)

From this general solution, the case of second order kinetics is obtained as a special case by setting ˇ = 0. The solution becomes now 2

b 1 + ˇcLA y

=0

y= √ 1 + 4˛ + 1

(43)

Again the dimensionless quantities are introduced: ˇ = b and ˛ =  (1 + ˇ ). After some straightforward algebraic ˇcLA manipulations, the physically meaningful analytical solution is obtained: y=

y=

2



2 (˛ + 1 − ˇ )

+ 4ˇ + ˛ + 1 − ˇ

(44)

1 ˛+1

(42)

(45)

Finally, the exceptional case of zero order kinetics is obtained by applying Eq. (39) into Eq. (29). The result becomes (46)

For the case of rapid mass transfer through the product layer, → 0, ˛ → 0 and y → 1 in Eqs. (44)–(46).

4.2.

Product layer model

1−y−

b )y  (1 + ˇcLA

y=1−

Rate equations (35) and (37) are obtained, for instance, from two-step reaction mechanisms, as has been shown in Murzin and Salmi (2005). The second and first order rate equations (36) and (38) are obtained as special cases of the more general Eqs. (37) and (35) by setting ˇ = 0.

4.1.

1−y−

For the case of first order kinetics, ˇ = 0 and the solution is obtained as a special case from Eq. (44):

The application of the proposed models on various rate equations, gives the surface concentrations. The following rate equations relevant for real fluid–solid processes will be considered: R =

For the limiting case that the mass transfer through the product layer is rapid compared to the chemical reaction rate, ˛ → 0 (since → 0) and y → 1 for both cases, Eqs. (41) and (42). For the rate equation (37), the following equation is obtained:

Shrinking particle model

For the rate equations (35)–(39) presented, the solution for the shrinking particle model can be obtained by setting  = 1 and replacing by  according to Eqs. (33) and (34). The parameter b . ˛ is now ˛ =  (1 + ˇ ), where ˇ = ˇcLA

4.3.

Summary of analytical solutions

The results for the product layer and shrinking particle models are summarized  s in Table 2. The values of the surface concentrations cLi needed in the solution of the local balance equations are obtained from Table 2 for the specific cases; for a general case, numerical solutions of Eqs. (28) and (32) are needed, at different reaction times and conversions of the solid substance. The modified equation helps especially in the reliable determination of the pre-exponential factor. The term Tmean represents the reference temperature of the experiments, k is the rate constant, T is the temperature (in Kelvin) and Ea is the activation energy.

5. Reaction area and reaction order with respect to the solid component For fluid–solid processes, the overall rate is assumed to be proportional to the accessible surface area (A). For instance, for a dissolved component (i) in a perfectly mixed, batchwise operated tank reactor, the balance equation can be written as

 dni ik Rk A = dt S

k=1

(51)

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

The area-to-volume ratio A0P /V0P is expressed by the shape factor (a) and the characteristic dimension (R0 ), a = (A0P /V0P )R0 , which gives

Table 2 – Simple kinetics: analytical solutions for product layer and shrinking particle models. 2 kcs LA 1+ˇcs LA kcs LA 1+ˇcs LA

R = R =

,

y= y=

R = k y = 1 − y=





2 (1−ˇ )2 +4(˛+ˇ )+1−ˇ 2

(˛+1−ˇ )2 +4ˇ +˛+1−ˇ

1883

(47)

A=

(48)

product layer and shrinking particle (49)

n0 M a AP P R0 A0P

(57)

Rearrangement of Eq. (57) gives

cs LA cb LA

aM AP n0 P R0 A0P

Rate constant: k = k0 e(−Ea /R)(1/T−1/Tmean ) (50)

A=

Explanations:

The ratio AP /A0P obtains the following values for the ideal geometries: (R/R0 )2 for a sphere, (R/R0 )1 for an infinitely long cylinder and (R/R0 )0 for a slab. After expressing these ratios with the particle volumes we obtain (VP /V0P )2/3 , (VP /V0P )1/2 and (VP /V0P )0 for the ideal geometries. The relation can be generalized to (VP /V0P )(a−1)/a for an arbitrary geometry. Eqs. (53) and (54) give

Product layer model See Eq. (27)

 

, ˇ ˛

b = −vi R (1) R/(DeA cLA )

 (1 + ˇ ) 



Shrinking particle model 1 



=

=−

vi R (1) (KLA cb ) LA

(1 + ˇ )



The modified equation helps especially in the reliable determination of the pre-exponential factor. The term Tmean represents the reference temperature of the experiments, k is the rate constant, T is the temperature (in Kelvin) and Ea is the activation energy.

The specific rates (Rik ) are expressed as mol/(m2 s). The totally accessible surface area (A) is related to the number of particles (nP ) and the specific surface area (). The specific surface area is a useful concept in this connection, since it can be measured by gas adsorption from samples withdrawn from the reaction mixture. The relation between the available particle surface area and the specific surface area for a general particle geometry will be derived in this section. For the sake of simplicity, equal particle sizes are presumed, but the treatment can be easily extended to comprise particle size distributions, too. By assuming equal-size particles, the accessible total surface area is A = nP AP

(58)

VP n = V0P n0

(59)

Consequently, we obtain AP = A0P

 V (a−1)/a P

V0P

(60)

and AP = A0P

 n (a−1)/a n0

(61)

which is inserted in Eq. (58): A=

aM n0 1/a n1−1/a P R0

(62)

The ratio a/R0 is a A0P = R0 V0P

(63)

(52) which can be expressed with the original specific surface area

The volume of a particle is = nM mP = VP =  nP P

A0P m0P

(64)

(53) Thus we have

where M and n denote the molar mass and the total molar amount, respectively. The number of particles can be obtained from the initial state (0),

a m0P = = P R0 V0P

(65)

which is inserted in the expression of A: V0P

n0 M = nP P

(54)

which gives nP =

n0 M V0P P

(55)

n0 M A0P AP P V0P A0P

(56)

(66)

If the mole fraction of the solid material is constant (x0j ), the amounts of substance are: n0j = x0j n0 (t = 0) and nj = x0j n (t > 0). Inserting these expressions in Eq. (66) gives A=

The accessible area becomes A=

A = Mn0 1/a n1−1/a

M x 1−x n n x0i 0i i

(67)

This expression for the area is used in the bulk-phase mass balances of the fluid and solid components.

1884

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

6. Mass balances for components in batch reactor

Table 3 – The effect of surface roughness and porosity on the observed reaction order.

Application of Eq. (67) to the mass balance of a solid component (j) gives (only one reaction is assumed to proceed in the system) dnj dt

=

j M x0j

n0j x nj 1−x R

(68)

R0 /␮m

a = 1/x

1−x

3 15 30 300

2/3 14/15 29/30 299/300

1 5 10 100

Note:  = 2.0 m2 /g,  = 1500 kg/m3 , a = R0 , x = 1/a, 1 − x = observed reaction order.

For a liquid-phase component (i) we obtain dni  M x 1−x  n n R = i dt x0j 0j j

(69)

It is worth noting that in spite of the complexity of the system, the very simple stoichiometric relationship is valid. Division of Eq. (68) by Eq. (69) gives the ratio of stoichiometric coefficients for the case that only one reaction proceeds in the system. Very often the component concentrations are used in the further treatment of the mass balances. For the liquid-phase components, the use of concentrations is evident. The amount of substance (ni ) is divided by the fluid volume. Formally a similar definition of the concentration can be applied to the solid component, too. The next step is crucial—is the fluid volume constant or not. For a varying fluid-phase volume (due to evaporation or chemical changes in the system), it is better to remain with the original forms of the balance equations (68) and (69). If the fluid volume can be considered to remain constant, the balance equations (68) and (69) become (70) and (71) after division by the fluid volume, dcj dt

=

j M x0j

x 1−x  c0j cj R

dci  M c0 cj 1−x R = i dt x0j

(70)

The concentrations at the reaction surface are calculated iteratively or analytically from Eqs. (28) and (32) for the product layer and shrinking particle models, respectively. In case of chemical rate control (rapid mass transfer), the bulk-phase b concentrations cLi are used directly in the rate expressions. A stiff ODE solver can be used to obtain the bulk-phase concentration profiles (Hindmarsh et al., 1980). In case of fitting the model to experimental data, nonlinear regression analysis is applied. The ordinary differential equations (ODE) describing the model are solved several times during the search of the best parameter values. The estimation of kinetic parameters can carried out by nonlinear regression analysis using the simulation and parameter estimation software MODEST (Haario, 2001). The system of ordinary differential equations is solved with the backward difference method implemented in the software ODESSA (Leis and Kramer, 1988), which is based on the stiff ODE solver, LSODE (Hindmarsh et al., 1980) code. The objective function, i.e. the sum of residual squares, Q, is minimized with the hybrid Simplex–Levenberg–Marquardt method,



2

Q = cexp − cest  =

nsets nobs(k)   nydata(j,k)  k=1

(71)

Typical rate expressions (R ) are shown in Table 2 (i = A). Eqs. (70) and (71) reveal an interesting concentration dependence of the observed overall rate (dcj /dt) with respect to the solid component. The term is constant throughout the reaction, while the concentration ci varies. The experimentally observed reaction order thus depends on the particle geometry. For a long cylinder, x = 1/2 and for a perfect, non-porous sphere x = 1/3. Thus reaction orders (=1 − x) 1/2 and 2/3 are predicted for these particles. For a slab, x = 1 and the observed reaction order approaches zero. However, as the real particles deviate from the ideal ones, the accessible surface area for the chemical reactions is higher. Parameter x has in general the value x = 1/a, in which a  3 in many cases. This implies that effective reaction orders 2/3 < 1 − x < 1 appear in practise. For very porous particles, x is very small and first order kinetics is observed experimentally. Thus, reaction orders exceeding 2/3 are in general very probable for real particles with some surface roughness. This issue is illustrated in Table 3. The table indicates that particles with larger radii than 5–10 ␮m easily exhibit first order behaviour with respect to the solid substance. In a general case, the concentration profiles as a function of the reaction time can be obtained by solving numerically the bulk-phase balance equations (68) and (69) or (70) and (71).

j=1

(cexp,ijk − cest,ijk )

2

(72)

i=1

where cexp represents experimental data (typically concentrations) and cest the estimated values (concentrations predicted by the model); nsets, nobs and nydata denote the number of data sets (number of batch experiments), number of experimental observations in the set and number of components, respectively. This implies that all experimentally recorded concentrations are incorporated in the minimization of the objective function.

7.

Numerical simulation example

In order to demonstrate how the generalized model works in practise, simulations were performed with values corresponding to commonly observed experimental data. The general model is able to account for varying morphology and composition of the solid phase, external- and internal mass transfer limitations as well as film diffusion resistance, variations in concentrations and temperatures, stoichiometry and varying reaction mechanisms and kinetics. The focus will be put on demonstrating the most important phenomena, as it is not practically possible to demonstrate all of the special cases in the current work. In the simulations, a solid phase is let to react with a liquid phase in different experimental conditions and one parameter is varied at a time to facilitate the influence on the model prediction. The values of the parameters used in modelling are displayed in Table 4.

1885

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

Table 4 – Parameter values used in numerical simulations. Parameter k Bi Dei R (m) 

Value −10

−8

Parameter

2

3

n

Value

2

1.0 –1.0 (mol/m )(m /mol) /s 10–100 5.0−10 –5.0−12 m2 /s (10–1000)10−6 −1

 (m /g) M (kg/mol) B x0j t (s)

Parameter

Value

3

1 0.001 0.001 1 0–3600

c0i (mol/m ) a c0j (mol/m3 )

2000 3–15 1000

2000 1900 1800

k=1.0 10-10

mol/m3

1700 1600 1500

k=5.5 10-10

1400 1300 1200 1100 1000

k=1.0 10-9 0

500

1000

1500 2000 time (s)

2500

3000

3500

Fig. 5 – The surface concentration (– – –) and the bulk concentration (—) for different values of the rate constant (product layer model). The lower curve in the pairs of curves represents the surface concentration.

8.

Experimental example: Zinc leaching

The experimental example is taken from metallurgical industry. Zinc appears as sulphide in ores and the sulphide is treated further by leaching. In the direct leaching of zinc concentrate with ferric iron, zinc is dissolved and ferric iron is simultaneously reduced to ferrous iron in an acidic environment by the sulphide in the concentrate according to the following overall reaction stoichiometry: Fe2 (SO4 )3 (aq) + ZnS(s) → 2FeSO4 (aq) + ZnSO4 (aq) + S(s)

1000 900 800 700

mol/m3

The importance of the relation between the reaction rate and the internal mass transfer limitations can be demonstrated by beginning with slower chemical kinetics and gradually increasing the chemical reaction rate in a case where a product layer is formed on the surface of the particle. The parameter values for external and film diffusion limitations are adjusted so that they are not limiting factors. Fig. 5 displays how the relation between the surface concentration and the bulk concentration changes when the reaction rate is increased. A dramatic difference can be observed due to internal diffusion limitations. The reaction kinetics is chemically controlled at the lowest reaction rate (surface concentration is equal to the bulk concentration), while internal mass transfer dominates heavily when the rate constant increases e.g. due to temperature. The shape of the particle influences the overall reaction rate very much, particularly at high conversions. Fig. 6 demonstrates how the overall reaction rate varies when the shape factor (a) is increased, i.e. when the surface becomes more rough. For high values of the shape factor, the reaction order of the solid component approaches one, and small solid rests remain unreacted. The shrinking particle model is used in Fig. 6, which corresponds with negligible internal mass transfer resistance. Depending on the experimental conditions and reaction media, the particle radius can also heavily influence the dissolution rate of the solid phase. This can easily be observed in the model prediction by changing the radius, while having a lower value (10) for the Biot number, which corresponds with film diffusion resistance. The other parameter values are kept as previously. Fig. 7 depicts how the surface concentration of the particle and the bulk concentration changes as a function of time. The reaction rate increases substantially when smaller particles are used, even though the specific surface area is kept constant.

600 500 400 300 200 3

100 0

0

500

1000

1500

5

2000

time (s)

9 2500

15

3000

3500

4000

Fig. 6 – The change in the concentration of the solid phase as a function of time for particles of varying morphology (roughness of the particle) in the absence of external mass transfer limitations (shrinking particle model). The values for the shape factor (a) are displayed in the figure.

1886

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

2000

r=1 mm

1800

r=0.5 mm

1600

mol/m3

1400

r=0.01 mm

r=0.01 mm

r=0.05 mm

1200 1000

r=0.05 mm

800 600 400

r=0.5 mm r=1 mm

200 0

0

500

1000

1500

2000

2500

3000

3500

4000

time (s) Fig. 7 – The surface concentration (– – –) and the bulk concentration (—) for different particle sizes as a function of time, when the Biot number is fixed to 10. In current commercial operations, the produced ferrous iron is re-oxidized with oxygen in order to continue the leaching. However, if the aim is to reduce the ferric iron concentration, the sphalerite ((Zn,Fe)S) concentrate can successfully be used as a reducing agent and is then leached at the same time. This would be a stage in a new leaching process, where jarosite (KFe3+ 3 (OH)6 (SO4 )2 ) is precipitated without a need for a neutralization step and thus, the problem with losses of zinc is avoided (Fugleberg, 2002). Reactive particles are illustrated in Fig. 4. Kinetic studies have been published for the abovementioned reaction, but the rate models proposed are often of empirical character and thus of very limited use. Competing hypotheses exist, whether the reaction proceeds on the solid surface according to the product layer model (Rath et al., 1981; Palencia Perez and Dutrizac, 1991; Kolodziej and Adamski, 1990; Suni et al., 1989; Crundwell, 1987) or in the liquid phase (Verbaan and Crundwell, 1986; Lotens and Wesker, 1987; Crundwell and Verbaan, 1987; Kammel et al., 1987). It is generally not agreed, if the shrinking particle or product layer model should be applied for zinc sulphide. Previous investigations do not strongly support the product layer (Markus et al., 2004), but more probably a shrinking particle model. Recently published research inclines more to the shrinking particle model (Xie et al., 2007; Dutrizac, 2006). Moreover, chemically (Dutrizac, 2006; Souza et al., 2007; Gökhan, 2009) and electrochemically (Verbaan and Crundwell, 1986) controlled surface reactions have been discussed in the literature and used for the modelling of the dissolution reaction. However, for pure ZnS as well as for sphalerite a chemically controlled surface reaction has been successfully applied to describe the leaching process in many studies. The contradictions in the mechanisms could be explained by the fact that models assuming spherical particles were used in several studies, even though the particles are very rough in morphology. Furthermore,

non-integer exponents for the ferric ion concentrations have been obtained to in the model fitting, which complicates the model discrimination. Several kinetic models were first tested on data obtained from a well-stirred isothermal laboratory-scale reactor. The models were obtained from Table 1. Soon it became clear that the classical shrinking particle and product layer models based on perfect, non-porous spheres could not describe the data well: They predicted the wrong reaction order with respect to the solid species. An empirical approach revealed that the reaction order with respect to the solid (Zn) species is very close to one, wile the reaction order with respect to the liquid species (Fe) exceeds one (it was around 1.25) (Grénman et al., 2007). In order to obtain a deeper understanding of the dissolution process, more rigorous modelling and evaluation of the obtained results was performed. For the reaction between ZnS(s) and ferric ions (Fe(III)), a stepwise surface reaction mechanism was proposed. The ferric ion forms a surface complex (I1 ) with zinc sulphide, the complex reacts further with another ferric ion forming a surface intermediate (I2 ) which releases ferrous ions (Fe(II)) and elementary sulphur is formed. Stepwise desorption of ferrous ions is also examined. The shrinking particle model was applied for zinc sulphide: Elementary sulphur is removed from the surface. Based on these considerations, the reaction mechanism can be written as ZnS(s) + Fe3+ ↔ I1

I1 + Fe3+ ↔ I2

I2 ↔ S(s) + 2Fe2+ + Zn2+ ZnS(s) + 2Fe3+ ↔ S(s) + 2Fe2+ + Zn2+

(I)

(II)

(III)

1887

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

0.2

75°C 85°C

Fe 3+ (mol/L)

0.15

95°C

0.1

0.05

0 20

0

40

60 Time (min)

80

100

120

Fig. 8 – Leaching of zinc with ferric iron: the fit of the model to experimental data at different temperatures; rate controlled by intrinsic kinetics. The solid lines represent model predictions. 0.2 200 rpm 350 rpm 500 rpm

Fe3+ (mol/L)

0.15

Application of the quasi-steady state hypothesis (Murzin and Salmi, 2005) on the intermediates (I1 and I2 ) and elimination of their concentrations from the rate equations gave finally the following rate expression r=

0.1

0.05

0 0

20

40

60 Time (min)

80

100

120

Fig. 9 – The fit of the model to experimental data at different agitation rates for experiments performed at 85 ◦ C. The solid lines represent model predictions.

2 2 c − cFeII k(cFeIII ZnII /K)

D

(73)

where D = k−1 k−2 + k−1 k+3 + k+2 k+3 cFeIII , k = k+1 k+2 k+3 , where k+1 = k+10 c* ZnS and K = k+1 k+2 k+3 /(k−1 k−2 k−3 ) k−3 = k−30 c* S . The products of the forward (k+ ) and backward (k− ) rate constants can be merged to single parameters and the overall reaction turned out to be irreversible (K has a high value). Thus the final, simplified form of the rate equation becomes 2

R =

s kcFeIII

(74)

s 1 + ˇcFeIII

Selection of the experimental system and equipment

Kinetic investigations

Structural investigations

Mass- and heat transfer studies

Ideas on the reaction mechanism including structural changes of the solid

Derivations (and simplification) of rate equations

Estimation of kinetic and mass transfer parameters

Model verification by numerical simulations and additional experiments

Fig. 10 – Scheme for investigation of fluid–solid reactions.

1888

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

which is one of the cases presented in Table 2; thus the surface concentration was calculated during the parameter estimation from the corresponding y-value (see Table 2): s cFeIII =



2 2

(1 − ˇ ) + 4(˛ + ˇ ) + 1 − ˇ

cFeIII

(75)

The morphology of the solid phase was implemented via the shape factor and a kinetic model was obtained for the process. The mass balance equations become dcZnS ZnS M x c c1−x R = dt x0ZnS 0ZnS ZnS

(76)

dcFeIII  M x c c1−x R = FeIIIi dt x0ZnS 0ZnS ZnS

(77)

The algorithm works as follows: the surface concentration is calculated from Eq. (75) and inserted into Eq. (74) to get R , which is inserted in Eqs. (76) and (77) to obtain the derivatives needed in the numerical solution of the ODEs. An activation energy of 53.2 kJ/mol was obtained for the dissolution. A good fit of the model to the experimental data was obtained in the absence and presence of mass transfer limitations. Examples of the model fits are displayed in Figs. 8 and 9. The present case exemplifies how the modelling procedure can help to determine the solid phase exponent correctly, which can further be verified by SEM. This leads to a much better result than forcing a model with an ideal morphology (e.g. sphere) to account for the data. Moreover, it illustrates how a reaction mechanism (in this case for a stepwise reaction) can be revealed by further modelling once a reliable exponential value is determined for the solid phase.

9.

Conclusions

The essential features of fluid–solid reaction systems were reviewed: reaction kinetics, internal and external mass transfer as well as surface morphology. New, generalized models for shrinking particles and product layers were proposed. The models are based on the use of a shape factor (a), which can obtain high and non-integer values (Table 3) for real particles with surface defects. The generalized models give a direct relation between the accessible surface area, shape factor and specific surface area (Eq. (67)). The modelling approach was illustrated with numerical simulations and an experimental example, leaching of zinc ores with ferric ions. Investigation of fluid–solid reactions requires an interdisciplinary approach, combining the studies of mechanism, kinetics, transport phenomena and particle morphology. The essential feature is the combination of experiments and modelling. The approach is shortly summarized in the scheme (Fig. 10).

Acknowledgement The work is a part of the activities of Process Chemistry Centre (PCC), a Centre of Excellence financed by Åbo Akademi. Financial support from Academy of Finland is gratefully acknowledged.

References Andersson, N., Ph.D. Thesis 2003. Modelling of Kraft Cooking Kinetics Using Near Infrared Spectroscopy. Karlstad University Studies, 2003:21, ISBN: 91-85019-36-4.

Borgen, J., Brelid, H., Theliander, H., 2007. Reaction kineticas of softwood kraft delignification—general considerations and experimental data. Nord. Pulp Pap. Res. J. 22, 177. Borgen, J., Brelid, H., Theliander, H., 2009. Towards a general kraft delignification model. Nord. Pulp Pap. Res. J. 24, 33. Christensen, T., Albright, L.F., Williams, T.J., 1983. A kinetic mathematical model for the kraft pulping of wood. TAPPI Proc., 239. Crundwell, F., 1987. Kinetics and mechanism of the oxidative dissolution of a zinc sulphide concentrate in ferric sulphate solution. Hydrometallurgy 19, 227. Crundwell, F., Verbaan, B., 1987. Kinetics and mechanism of the non-oxidative dissolution of sphalerite (zinc sulphide). Hydrometallurgy 17, 369. Debasish, S., Bandyopadhyay, A., 2011. Shrinking core model in characterizing aqueous phase dye adsorption. Chem. Eng. Res. Des 89, 69. Dickinsin, C., Heal, G., 1999. Solid–liquid diffusion controlled rate equations. Thermochim. Acta 340-341, 89. Dutrizac, J.E., 2006. The dissolution of sphalerite in ferric sulfate media. Metall. Mater. Trans. B 37B, 161. Fabiano, M., Pina, P., Porcaro, R., Oliveira, V., Silva, C., Leao, V., 2010. The kinetics of zinc silicate leaching in sodium hydroxide. Hydrometallurgy 102, 43. Fugleberg, S., Patent No. WO0246481 2002. Method for the Hydrolytic Precipitation of Iron. Gökhan, U., 2009. Kinetics of sphalerite dissolution by sodium chlorate in hydrochloric acid. Hydrometallurgy 95, 39. Grénman, H., Murzina, E., Rönnholm, M., Eränen, K., Mikkola, J.-P., Lahtinen, M., Salmi, T., Yu, Murzin D., 2007. Enhancement of solid dissolution by ultrasound. Chem. Eng. Process. 46 (9), 862–869. Grénman, H., Eränen, K., Krogell, J., Willför, S., Salmi, T., Murzin D.Yu., 2011. The kinetics of aqueous extraction of hemicelluloses from spruce in an intensified reactor system. Ind. Eng. Chem. Res. 50, 3818–3828. Grénman, H., Wärnå, J., Mikkola, J.-P., Sifontes, V., Fardim, P., Murzin D.Yu. Salmi, T., 2010. Modeling the influence of wood anisotropy and internal diffusion on delignification kinetics. Ind. Eng. Chem. Res. 49, 9703–9711. Haario, H., 2001. Modest 6.0-User’s Guide. Profmath Oy, Helsinki. Habashi, F., 2005. A short history of hydrometallurgy. Hydrometallurgy 79, 15. Hindmarsh, Lsode, A., Lsodi, 1980. Two new initial value ordinary differential equation solvers. ACM SIGNUM Newslett. 15, 10. Hossein, A., Fazlollahi, F., Teherani, S., 2010. Leaching kinetics of calcined magnesite in ammonium chloride solution. Aust. J. Basic Appl. Sci. 4, 5956. Kammel, R., Pawlek, F., Simon, M., Xi-Ming, L., 1987. Oxidizing leaching of sphalerite under atmospheric pressure. Metall 41, 158. Kolodziej, B., Adamski, Z., 1990. Dissolution of sphalerite in aqueous hydrochloric acid solutions under reduction conditions. Hydrometallurgy 24, 393. Leis, J., Kramer, M., 1988. Odessa—an ordinary differential equation solver with explicit simultaneous sensitivity analysis. ACM TOMS 14, 66. Levenspiel, O., 1999. Chemical Reaction Engineering, third ed. Wiley, New York. Lotens, J., Wesker, E., 1987. The behavior of sulfur in the oxidative leaching of sulfidic minerals. Hydrometallurgy 18, 39. Markus, H., Fugleberg, S., Valtakari, D., Salmi, T., Murzin, D., Yu Lahtinen, M., 2004. Reduction of ferric to ferrous with sphalerite concentrate, kinetic modelling. Hydrometallurgy 73, 269. Murzin, D., Salmi, T., 2005. Catalytic Kinetics. Elsevier, Amsterdam. Órfão, J., Martins, F., 2002. Kinetic analysis of thermogravimetric data obtained under linear temperature programming—a method based on calculations of the temperature integral by interpolation. Thermochim. Acta 390, 195.

chemical engineering research and design 9 1 ( 2 0 1 3 ) 1876–1889

Palencia Perez, I., Dutrizac, J., 1991. The effect of the iron content of sphalerite on its rate of dissolution in ferric sulphate and ferric chloride media. Hydrometallurgy 26, 211. Rath, P., Paramguru, R., Jena, P., 1981. Kinetics of dissolution of zinc sulphide in aqueous ferric chloride solution. Hydrometallurgy 6, 219. RISI Report, 2009. World Market Pulp Capacity Report—Excerpt. Santos, A., Rodrigues, F., Gilarranz, M., Moreno, D., Garcia-Ochoa, F., 1997. Kinetic modeling of kraft delignification of eucalyptus globulus. Ind. Eng. Chem. Res. 36, 4114. Souza, A.D., Pina, P.S., Leao, V.A., Silva, C.A., Siqueira, P.F., 2007. The leaching kinetics of zinc sulphide concentrate in acid ferric sulphate. Hydrometallurgy 89, 72. Suchun, Z., Nicol, M., 2010. Kinetics of dissolution of ilmenite in sulphuric acid solutions under reducing conditions. Hydrometallurgy 103, 196. Suni, J., Henein, H., Warren, G., Reddy, D., 1989. Modelling the leaching kinetics of a sphalerite concentrate size distribution in ferric chloride solution. Hydrometallurgy 22, 25.

1889

Szekely, J., Evans, J.W., Sohn, H.Y., 1976. Gas–Solid Reactions. Academic Press, New York. Trambouze, P., van Landeghem, H., Wauquier, J.P., 1988. Chemical Reactors—Design/Engineering/Operation. Editors Technip, Paris. Verbaan, B., Crundwell, F., 1986. An electrochemical model for the leaching of sphalerite concentrate. Hydrometallurgy 16, 345. Willför, Holmbom, S., Isolation, B., 2004. Characterization of water-soluble polysaccharides from Norway spruce and Scots pine. Wood Sci. Technol. 38, 173. Willför, S., Sundberg, A., Hemming, J., Holmbom, B., 2005. Polysaccharides in some industrially important softwood species. Wood Sci. Technol. 39, 245. Xie, K-q., Yang, X-w., Wang, J-k., Yan, J-f., Shen, Q.-f., 2007. Kinetic study on pressure leaching of high Iron sphalerite concentrate. Trans. Nonferrous Metal Soc. China 17, 187. Yang, L., Liu, S., 2005. Kinetic model for kraft pulping process. Ind. Eng. Chem. Res. 44, 7078.