A Versatile Lattice Simulator for FluidâSolid Noncatalytic Reactions in

Ind. Eng. Chem. Res. 1997, 36, 4993-5009

4993

A Versatile Lattice Simulator for Fluid-Solid Noncatalytic Reactions in Complex Media Alessandra Adrover*,†,§ and Massimiliano Giona‡,§ Dipartimento di Ingegneria Chimica, Universita` di Roma “La Sapienza”, Via Eudossiana 18, 00184 Roma, Italy, Dipartimento di Ingegneria Chimica, Universita` di Cagliari, and Piazza d’Armi, 09123 Cagliari, Italy, and Centro Interuniversitario sui Sistemi Disordinati e sui Frattali nell’Ingegneria Chimica, Universita´ di Roma, Via Eudossiana 18, 00184 Roma, Italy

A versatile lattice simulator for fluid-solid noncatalytic reactions is developed in detail in order to study linear, nonlinear, and nonisothermal kinetics in complex media, in the presence of multicomponent diffusion and n-solid reactant species. The simulator is based on a time-splitting algorithm for the diffusion and for the reaction steps. The simulator is quantitatively checked in many different cases involving initially nonporous particles constituted by a solid reactant dispersed in an inert matrix. The influence of spatial correlation properties in the case of uniform and nonuniform solid reactant distribution is analyzed in detail. 1. Introduction Fluid-solid reactions represent a basic mechanism in many industrial processes, from coal gasification to ore leaching. Over the last 20 years, many authors have studied the application of percolation models and of other concepts derived from the theory of disordered systems to interpret the relationship between the macroscopic overall kinetics and structural properties of the solid matrices, especially in regard to porous solids (Sahimi, 1994). The first attempt to apply percolation models was performed by Mohanty et al. (1982). Subsequently, Reyes and Jensen (1986a,b, 1987) modeled different gas-solid reaction (char gasification, sulfation of calcined limestones) in both reaction-controlled and diffusion-controlled regimes. The estimate of the structural properties and of effective transport parameters comes from the assumption of a Cayley-tree topology in the pore network (Yortsos and Sharma, 1986; Kerstein and Bug, 1986). These models are effective continuum models, in which the balance equations for the fluid reactant are written for continua, and the effective transport parameters and the macroscopic description (and temporal evolution) of the pore network (accessible porosity, accessible surface area) are obtained from percolation models or from other models of disorder. There are several drawbacks to continuum models for disordered media. The first is that they may not be valid in the presence of structures possessing infinitecorrelation length. This occurs, for example, in gassolid reactions in the proximity of the fragmentation threshold (since fragmentation can be regarded as a critical phenomenon). Moreover, all the parameters (transport and structural) entering into hybrid models are based on a priori assumptions about the pore-network structure and/or distribution of the solid reactant within an inert matrix. In particular, transport and structural parameters are * Author to whom correspondence should be addressed. Phone: +39-6-44585892. Fax: +39-6-44585339. E-mail: alex@ giona2.ing.uniroma1.it. † Universita ` di Roma “La Sapienza”. ‡ Universita ` di Cagliari. § Centro Interuniveristario sui Sistemi Disordinati e sui Frattali nell’Ingegneria Chimica. S0888-5885(97)00164-4 CCC: $14.00

generally obtained starting from simple (and uncorrelated) models of disorder (the classical site and/or bond percolation models) or from unrealistic topologies (such as Cayley trees, which are loopless and nonfinite in dimension; see Aharony and Stauffer (1992)). The limitation of these assumptions is that both the porenetwork structure of real materials and/or the distribution of solid reactants within the pore-network matrix are spatially correlated in real materials. Although the spatial correlations are usually short-range (and the universality class is therefore the same as the corresponding uncorrelated structures), the structural properties, for noncritical structures (i.e., not at the percolation threshold), may be different, as shown by Adrover and Giona (1996) for permeability. A detailed discussion of this issue is developed in Section 5. Over the last 10 years, the growth of interest in the influence of the structural features of porous and granular materials on global transport and reaction properties together with the availability of fast computers possessing large memory capacity has stimulated the development of lattice models of porous structures. Lattice models for noncatalytic fluid-solid reactions were developed by Kerstein and Edwards (1987) and by Sahimi and Tsotsis (1987, 1988). Kerstein and Edwards (1987) modeled the char oxidation process as the ignition and burnout (i.e., removal) of solid bonds. Upon ignition, a solid bond is assigned a randomly selected burning time drawn from a distribution function depending on the mass of the fragment. The distinguishing features of the film-diffusion-limited regime of char oxidation are incorporated into the solidbond ignition criterion and into the dependence of the particle-burning rate per unit surface area on the mass (i.e., the number of solid bonds) of a given subcluster. In the model developed by Sahimi and Tsotsis (1987, 1988), diffusion is simulated by means of a random walk of the fluid reactant, which is controlled by the local bond conductances (diffusivities), and reaction is a random event occurring with a probability depending on the rate of reaction (and on temperature). The simulation refers, however, to isothermal conditions. This kind of model takes advantage of the basic analogy between transport (diffusion) and random walk, as widely applied in the literature on elementary reactions in disordered (fractal) media. Lattice models can be classified in terms of two major © 1997 American Chemical Society

4994 Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997

categories: effective lattice models and strictly geometrical models. The basic prototype for effective lattice models is given by the dual site-bond model proposed by Mayagoitia and co-workers (Mayagoitia et al., 1989a; Javier-Cruz et al., 1989) for porous structures and extended by Zgrablich and co-workers (Riccardo et al., 1988, 1992; Mayagoitia et al., 1989b, 1990) to cover the case of energetically heterogeneous adsorbents (Zgrablich et al, 1996; for a review see Zgrablich, 1997 and references therein). The structure of a porous medium is regarded as a network of sites, corresponding to the voids of the pore space, and of bonds, corresponding to the necks connecting the sites. Each solid bond is assigned a geometric radius in accordance with two main principles, as discussed in Mayagoitia et al. (1989a) and Javier-Cruz et al. (1989). The parametrization with respect to the radii is particularly suitable to describe physical phenomenologies, such as N2-adsorption hysteresis, mercury porosimetry, two-phase immiscible displacement, etc. (Faccio et al., 1993; Sapag et al., 1993; Gonzales et al., 1994), the physics of which is controlled by capillarity and wettability phenomena described by such equations as the Kelvin equation or the Laplace equation, which depends on the local curvature (pore radius). These models do not correspond to a specific geometric structure, i.e., to a set of points in the space occupied either by the pore network or by the solid matrix. The other class consists of strictly geometrical models in which a particle is viewed as a set of points described in the most complete (pointwise) way by means of a characteristic function χP (x) which, for each lattice site x, returns either 0, if site x belongs to the pore space, or an integer i ) 1, ..., N, if x corresponds to the ith species in the solid mixture. The prototypes of any geometric lattice model are sitepercolation clusters (Aharony and Stauffer, 1992). More realistic descriptions of real porous materials have been developed by Joshi (1974), Quiblier (1984), Adler et al. (1990), and Giona and Adrover (1996) to account for the effects of short-range spatial correlations. Within the framework of transport and reaction phenomena, strictly geometrical lattice models are particularly appealing due to the local nature of the reactive events and to nearest-neighboring interactions in the description of transport. However, the lattice models developed so far suffer from the intrinsic limitation that they are applicable only to elementary reactions and to simple transport models (pure diffusion with constant diffusivities). The great merit of the lattice models, based on random walk simulation and on simple rule to implement surface reactions, is that they envisage the qualitative features and the anomalies associated with the interplay between network topology and elementary transport/reaction kinetics (Havlin and Ben-Avraham, 1987). Nevertheless, none of these models can really describe more complex kinetics, such as nonisothermal kinetics in the presence of simultaneous heat and mass transfer, or reactions in the presence of multicomponent transport phenomena in the fluid phase. These limitations of random-walk simulations have been stressed by Giona and Viola (1994) and Giona et al. (1994) in the study of nonlinear diffusion equations. In the case of concentration-dependent Fickian diffusion coefficients, which occur in the overwhelming majority of multicomponent transport problems (Taylor and

Krishna, 1993), it is of course still possible to formulate Monte Carlo simulation schemes. As a result of nonlinearity, however, Monte Carlo methods prove to be quite inefficient since they would involve the simultaneous analysis of a statistically significant ensemble of random walkers in order to account for the concentration-dependent effects. If we are to achieve a complete understanding of the relationship between structural properties and overall kinetic behavior in a way that could be useful for engineering applications of disordered models to chemical reaction theory, it is important to make use of a geometric lattice characterization of the solid matrix and to achieve an efficient and general way of representing the reactive and diffusion steps for arbitrary kinetics and transport conditions. One way to bridge the gap between approximate continuum theories and simplified random-walk simulations is provided by versatile lattice simulators, which are capable of overcoming the limitations of the network models so far developed. Such a versatile simulator is required to fulfill several basic requirements: (1) It should be applicable to generic lattice structures constituted by a given number of solid reactants, an inert solid matrix, and pores and to more than one fluid reactant. (2) It should give consistent results in both reaction-controlled and diffusion-controlled conditions, in the case of both initially nonporous particles and porous particles. (3) It should be applicable to arbitrary intrinsic kinetics: linear, nonlinear, and nonisothermal. (4) In those cases in which macroscopic continuum models can be applied, (e.g., the case of initially nonporous particles by considering the shrinking core (SC) model with the effective transport parameters characterizing the system), the simulator should be quantitatively in agreement with the continuum solutions over the entire range of operating conditions. In this article we propose a lattice simulator fulfilling all the requirements stated above. The simulator is based on a time-splitting algorithm (for diffusive dynamics and for surface reaction) and can be applied to generic lattice structures of arbitrary complexity. The diffusive dynamics is described by means of a simple and efficient algorithm encompassing zero-flux conditions at each solid boundary site (Giona et al., 1994, 1995). A brief explanation of the algorithm in a numerically extreme situation (i.e., for a fractal pore space) is given in the Appendix. Great care is used in testing the accuracy and the validity of the simulator in all cases of practical interest. Throughout the article, initially nonporous particles are considered, although the simulator can be applied with no change to porous solids (as discussed in Section 4.3). The reason for this choice is the following. In the case of initially nonporous structures, the Shrinking Core model can be applied in the case of moderate disorder (i.e., far from criticality) and in the diffusion-controlled regime (the most interesting case from the algorithmic point of view). Under these conditions, simulation results can be checked against the prediction of continuum models. Since the aim of this article is to envisage a lattice method to study reactions on disordered model structures, to be used in cases of practical interest, reactions on initially nonporous particles are the best benchmark for quantitative validation and for comparison with generally accepted results. From the computational point of view, the case of

Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997 4995

initially porous solids adds nothing new, while in a phenomenological framework, it is so rich and structured (since reaction evolution depends on the topology of the pore network, on the distribution of pore sizes and reactant solid clusters, etc.) that this case alone would necessarily require a detailed analysis in itself. This article is organized as follows. In Section 2 we develop the basic principles of the lattice simulator for generic nonlinear reactions. Section 3 discusses a manifold of different phenomenologies of increasing complexity and compares the simulation results with macroscopic continuum models. Section 4 extends the method to the general case of several solid and fluid reactants and to porous matrices. Section 5 analyzes a physical case study, namely the effects of spatial correlations and of radial disuniformities in reactant distribution, in order to show the effects of spatial correlation properties on the overall macroscopic kinetics. 2. The Algorithm The class of fluid-solid reactions considered throughout the article can be expressed by the following relation:

A(f) + bB(s) f C(f) + dD(s)

(1)

The pellet is initially nonporous and consists of a reactant and inert solid, so that the overall pellet size L and shape do not change with reaction (Sohn and Xia, 1986). The solid product layer is porous, the fluid species A diffuses through it with diffusivity DAC, and simultaneously a chemical reaction occurs on the unreacted solid boundary, with intrinsic isothermal kinetics RA [mol/s m2]:

-RA(cA) ) k(cA0)cAψ(cA/cA0)

(2)

where cA is the molar concentration of A, cA0 the external uniform fluid concentration, and ψ(cA/cA0) a dimensionless functional accounting for nonlinearities in the intrinsic kinetics. The concentration cA obeys the local diffusion equation

∂cA ) ∇(DAC∇cA) ∂t

(3)

equipped with the boundary conditions

n‚(-DAC∇cA) )

{

-RA(cA) on SR 0 on SI

(4)

SR and SI being respectively the solid reactant and the inert matrix surfaces and n a unit vector normal to the surface. Due to the reaction, the solid reactant boundary evolves with time. Let x ) x(t) be any point lying on the solid reactant surface in contact with the fluid species A and w(x,t) ) n‚x the local displacement in the direction orthogonal to the fluid-solid interface (Bekri et al., 1995). The displacement rate ∂w/∂t is related to the intrinsic kinetics RA through the relation

FB

∂w(x,t) ) -bk(cA0)cAψ(cA/cA0), ∂t

x ∈ SR

(5)

where FB is the molar density of the pure reactant B and b is the stoichiometric coefficient (eq 1). By introducing the dimensionless variables

η ) x/L, ξ ) cA/cA0, τ ) tk(cA0)/L

(6)

in the case of negligible external mass-transfer resistance (cA0 ) cA(x ) L)), the dimensionless characteristic parameters controlling reaction evolution are the square surface Thiele modulus φ2 ) k(cA0)L/DAC and R ) bcA0/ F B. The underlying principles involved in the construction of the lattice simulator are (1) the discrete description of the lattice geometry (solid/pore matrix) and (2) the continuous description of the concentration fields of solid and fluid reactants. This double level in the description of the structure enables us to apply the simulator to arbitrarily disordered lattice structures but also to avoid (or rather to control) the approximations typical of fully discrete lattice simulations, such as those where a solid reactant site disappears at each reactive event (i.e., the normalized solid concentration is assumed to be either 0 or 1). In a two-dimensional lattice representation of the pellet (in slab or cylindrical geometry), a characteristic function χij is defined for each lattice site (i, j), such that

χij )

{

1 if (i, j) ∈ P 0 if (i, j) ∈ B ∪ C

(7)

where P, B, and C are respectively the sets of sites representing the product layer (through which the fluid species diffuses), the reactant solid matrix, and the inert matrix (both impermeable to mass transport) [in a similar way, for a three-dimensional representation of the pellet, the characteristic function χijk is defined for each lattice site (i, j, k), such that χijk ) 1 if (i, j, k) ∈ P and χijk ) 0 if (i, j, k) ∈ B ∪ C ]; see Figure 1A. A continuous concentration field ξij is defined at all the sites belonging to the product layer, representing the dimensionless concentration of the diffusing species A. Moreover, a continuous concentration field qij is defined at all the solid reactant sites, qij, representing the unreacted fraction of the solid reactant at site (i, j) ∈ B. The lattice length scale is ∆η ) 1/N, N being the lattice size expressed in lattice units. A fundamental question in setting up of the simulator is the proper choice of the lattice time scale ∆τ′. Let ∆t be the time interval required for the conversion of ∆nB ) FB∆xd moles of solid reactant at the unit lattice site (d ) 2 for slab or cylindrical geometry, d ) 3 for spherical pellets). By taking into account the stoichiometry and the intrinsic isothermal kinetics, eq 2, it follows that [equation 8 can be seen as the discrete version of eq 5]

∆nB ∆xd ) FB ) bk(cA0)cAψ(cA/cA0)∆xd-1 ∆t ∆t

(8)

In dimensionless form, eq 8 implies that

∆η ) Rψ(ξ)ξ ∆τ

(9)

We may therefore define the lattice time scale ∆τ′ as

∆τ′ )

∆η Rψ(1)

(10)

which corresponds to the dimensionless time interval required for the advancement of the local reaction front at the maximun fluid concentration ξ ) 1 by the lattice


integer exceeding x) and ω ) ∆τs/(φ2∆η2Nd) < 1. Of course, ξij(0) ) ξij(τ) and ξij(Nd) ) ξij(τ + ∆τs/2). The finite-difference algorithm presented above can easily be extended to the case of position-dependent or concentration-dependent diffusion coefficients. In the case of a generic concentration-dependent diffusion coefficient DAC(ξ) ) DACof(ξ), eq 12 can be reformulated as

ξij(n+1) ) ξij(n) + ω 4f Figure 1. (left) Schematic representation of a slab pellet. (right) Pictorial representation of the Zij geometric factor for a cylindrical pellet.

length scale ∆η. For a slab geometry of the pellet, ∆τ′ coincides with the dimensionless time interval ∆τ′′ required for the complete conversion of a solid reactant site in contact with the fluid species A at maximum concentration ξ ) 1. Moreover, since the lattice simulator introduces a double level of description of the structure (a lattice representation of the solid reactant matrix and a continuous concentration field qij at all the solid reactant sites), it is possible to define the sampling time interval ∆τs as the fraction f < 1 of the time interval, i.e., ∆τs ) f∆τ. The utility and the choice of the parameter f are discussed below. During each sampling time interval ∆τs, the evolution of the sets P and B, and of the concentration fields ξij and qij, develops in two steps following a time-splitting procedure (Oran and Boris, 1986). The Diffusion Step. In the diffusion step, the concentration field ξij in the product layer is updated according to the diffusional dynamics of the fluid species A with the pellet. This step is essentially an explicit finite difference algorithm. In principle, the updating of ξij(τ) is given by

ξij(τ + ∆τs/2) ) ξij(τ) +

∆τs 2

φ ∆η

2

∑

4

∑

χhk[ξhk(n) - ξij(n)],

(h,k)∈I(i,j)

n ) 0, ..., Nd (12) where Nd ) [4∆τs/(∆η2φ2)] ([x] indicates the smallest

[ξhk(n) - ξij(n)] (13) χhk (n) (n) f(ξhk ) + f(ξij )

to be repeated Nd times, where Nd ) [4(∆τs/∆η2)f/φ2], ω ) 4(∆τs/∆η2)f/(φ2Nd), φ2 ) k(cA0)L/DACo, and f ) max{f(ξ)}, 0 e ξ e 1. The analogous of eq 11 in the case of concentration-dependent diffusion coefficient is discussed in the Appendix. The Reacting Step. At the beginning of the reacting step the fluid reactant concentration field is {ξij(τ + ∆τs/ 2)}. In the reacting step each site (i,j) belonging to B is analyzed by scanning the entire lattice. Let (i,j) be a solid reactant site in contact with the fluid species A (i.e., such that there exists a site (h, k) of I(i, j) belonging to P). The average fluid concentration 〈ξij〉 in the neighborhood of (i,j) is defined as 〈ξij〉 ) (1/νij)Σ(h,k)∈I(i,j)χhkξhk, with νij ) Σ(h,k)∈I(i,j)χhk. Due to the definition of the characteristic time interval ∆τ′ (see eq 10), in cases where the solid reactant site (i,j) is in contact with an average fluid concentration 〈ξij〉, only a fraction f∆qij of it is converted in the sampling time interval ∆τs ) f∆τ′. By using ∆τij ) ∆η/[Rψ(〈ξij〉)〈ξij〉] to indicate the dimensionless time interval required for the advancement of the local reaction front at (i,j) by ∆η, in the presence of the average fluid concentration 〈ξij〉, the quantity ∆qij is given by

∆qij )

(h,k)∈I(i,j)

where I(i,j) is the set of nearest-neighboring sites of site (i, j). This algorithm has already been applied by Giona et al. (1994, 1995) to study diffusion in the presence of complex boundaries, directly implementing the zero-flux local boundary condition on a solid site (see Appendix). Indeed, if site (h,k) ∈ I(i,j) does not belong to P, then the local flux J(h,k)f(i,j) ∼ χhk(ξhk - ξij) is automatically set to 0 since χhk ) 0. It should be observed that the factor ∆τs/φ2∆η2 can easily exceed 1/2 (which is the upper boundary for the stability of explicit finite-difference algorithms; see Crank (1967)), since the length and time scales [indeed ∆τs depends on the additional parameter f, the optimal choice of which is discussed below] are already fixed. This problem can be overcome by iterating the stable algorithm Nd times:

ω

(h,k)∈I(i,j)

χhk[ξhk(τ) ξij(τ)] (11)

ξij(n+1) ) ξij(n) +

∑

2f(ξhk(n))f(ξij(n))

∆τ′ 〈ξij〉ψ(〈ξij〉) ) ∆τij ψ(1)

(14)

In terms of the concentration field qij, eq 14 implies that

qij(τ + ∆τs) ) qij(τ) - f∆qij

(15)

In order to consume f∆qijFB∆x2 mol of the solid reactant site (i, j), by the enforcement of the stoichometry, f∆qijFB∆x2/b mol of the fluid species A in the neighborhood of (i, j) are needed. In dimensionless term, this means that the average fluid concentration 〈ξij〉 in the neighborhood of the reactant site must be updated as follows 〈ξij〉 f 〈ξij〉 - f∆qij/R. This is obtained by updating the concentration ξhk of the fluid species A at all the sites (h,k) ∈ P adjacent to (i,j) in the following way:

(

)

f ∆qij ξhk f ξhk 1 - χhk R 〈ξij〉νij

(h,k) ∈ I(i,j) ∩ P

(16)

It should be observed that, for low values of the R parameter, the quantity ∆qij/[〈ξij〉νijR] could exceed 1, leading to unphysical negative concentrations of A; eq 16. As a matter of fact, f can be chosen so that f/R < 1, in order to ensure the stability of the reaction step. The


possibility of tuning the value of the parameter f between 0 and 1 enables us to treat the whole range of R ∈ (10-3, 1), corresponding to gas-solid systems (R = 10-3) and liquid-solid systems (R = 0.1). When f∆qij g qij(τ), the reacting solid site (i, j) is eliminated and replaced by a new product-layer site. This implies that the sets B and P, and correspondingly the characteristic function χij, evolve in time. When ∆qij > qij(τ), ∆qij ) qij(τ) is assumed in the updating of the fluid species concentration (eq 16), in order to enforce stoichiometric balance. Such a discretization effect is practically negligible for low values of the parameter f (f e 0.05). Due to eq 16, at the end of the reaction step, the final fluid concentration ξhk(τ + ∆τs) is given by

∏

ξhk(τ + ∆τs) ) ξhk(τ)

(i,j)∈I(h,k)

(

)

f ∆qij 1 - χhk R 〈ξij〉νij

(17)

The updating of the concentration field in the reaction step depends in principle on the order in which the solid sites are visited. However, this effect is in absolutely negligible all practical cases, especially for small values of f. For a cylindrical pellet, eq 14 should be slightly modified in order to allow for the fact that the total conversion of a solid reactant site (i, j) corresponds to the local advancement of the reaction front (with respect to the center of symmetry of the pellet) by ∆rij ) ∆η/Zij (see Figure 1B), where

{

Zij ) max

xi2

,

}

xj2

(18)

xi2 + j2 xi2 + j2

Let us use ∆τ′′ to indicate the dimensionless time interval required for the complete conversion of a solid reactant site (i, j) at the maximum fluid concentration ξ ) 1. For a slab pellet ∆τ′ ) ∆τ′′. For a cylindrical pellet ∆τ′ ) ∆τ′′Zij and depends on the position of the site (i, j) itself. For this reason, the lattice time scale ∆τ′ is still given by eq 10, but this geometric effect is taken into account locally by rescaling (on the basis of the angular factors Zij defined by eq 18) the fraction of the solid reactant site ∆qij converted, i.e., eq 14 should be replaced by

∆qij ) [〈ξij〉ψ(〈ξij〉)]/[ψ(1)Zij]

{

xi2

,

xj2

,

mal kinetics, on ordered and disordered structures (percolation lattices) in order to show the validity of the simulator discussed in Section 2. Initially nonporous particles are considered so that the numerical results can be quantitatively compared with the corresponding SC model for continua. Slab and cylindrical pellets are considered, although the three-dimensional simulations can be performed exactly in the same way. 3.1. Linear Intrinsic Kinetics. Quantitative verification of the applicability and the validity of the lattice simulator can be obtained by first addressing the case of a linear isothermal kinetics (ψ(ξ) ) 1) for an initially nonporous particle (slab geometry) formed by the pure solid reactant B. The fluid species A diffuses through the product layer with constant diffusivity DAC. The product layer progressively fills the space initially occupied by the solid reactant in such a way that the overall pellet size and shape do not change with the reaction. Figure 2 shows the behavior of the total conversion X vs τ obtained by applying the lattice simulator (continuous lines), for different values of the surface Thiele modulus φ, compared with the solution of the dynamic (SC) model (dots):

1 ∂ d-1 ∂ξ ∂ξ ) η ∂τ ηdφ˜ 2 ∂η ∂η

(

(19) ξ ) 1 at η ) 1,

and the concentration fields qij and ξij are still updated in accordance with eqs 15 and 16. Analogous considerations hold for spherical pellets. In this case,

Zijk ) max

Figure 2. X(τ) vs τ obtained by the lattice simulator (continuous lines, N ) 200, R ) 1.0) for different values of the Thiele modulus φ, compared with the solution of the dynamic shrinking-core (SC) model (dots) with φ˜ ) φ. (a) φ2 ) 0.1. (b) φ2 ) 1.0. (c) φ2 ) 4.0. Dotted line represents the PSS solution of the SC model for R ) 1.0 and φ2 ) 4, i.e., in the case of significant deviation from the dynamic solution.

xk2

}

xi2 + j2 + k2 xi2 + j2 + k2 xi2 + j2 + k2

(20) At the end of the reacting step, the simulator furnishes a pointwise description of the evolving solid reactant matrix. The value of the total conversion can thus be obtained as X(τ) ) 1 - (1/NB)Σ(i,j)∈Bqij(t), where NB is the number of solid reactant sites initially present. 3. Results and Validation of the Algorithm This section analyzes different examples of increasing complexity, from linear intrinsic kinetics to nonisother-

)

ηc e η e 1

(21)

∂ξ ) φ˜ 2ξψ(ξ) at η ) ηc (22) ∂η

dηc ) -Rξψ(ξ)|η)ηc dτ

(23)

d ) 1, 2, and 3 respectively for slab, cylindrical, and spherical pellets, which can be properly applied in this case (by setting ψ(ξ) ) 1 and φ˜ ) φ) in both the reactioncontrolled (φ < 0.1) regime and the diffusion-controlled (φ > 1) regime. The adjective dynamic here refers to the solution of the unsteady-state shrinking-core model obtained by applying the method envisaged by Selim and Seagrave (1973a,b). The excellent agreement between simulations and the continuum solution, eqs 21-23, proves the correctness of the choice of the lattice time scale ∆τ′ and the validity of the procedure implementing the reaction step. The dotted line appearing in Figure 2 is the pseudosteady-state (PSS) solution of the SC model for R ) 1.0


Figure 3. X(τ) vs τ obtained by the lattice simulator (continuous lines, N ) 200) for initially nonporous slab particles formed by the solid reactant B randomly dispersed in an inert matrix, for different values of the uniform solid reactant fraction , in the diffusion controlled regime (φ2 ) 4). The intrinsic kinetics is linear, and the R value (R ) 0.1) corresponds to a liquid-solid system. (a) ) 1.0. (b) ) 0.9. (c) ) 0.8. (d) ) 0.7. Dots represent the solution of the dynamic SC model with φ˜ 2 ) φ2/DΣ().

and φ2 ) 4, i.e., in cases where the PSS solution diverges significantly from the dynamic solution. The PSS result is shown in order to highlight the fact that the lattice simulation (and specifically the implementation of the diffusion step) is able to give a correct reproduction of the unsteady-state diffusive dynamics of the fluid species within the product layer. Let us now turn our attention to disordered systems by considering initially nonporous pellets constituted by the solid reactant B randomly distributed within an inert matrix. The solid reactant distribution is generated by means of random (uncorrelated) site-percolation models (Aharony and Stauffer, 1992) characterized by the probability p, where here it represents the probability that a site belongs to the reactant matrix. For values of p higher than the critical threshold pc ) 0.5927, such as those considered in this section (p g 0.70), p practically coincides with the fraction of sites initially belonging to the reactant matrix. The initial lattice structure of the solid reactant matrix is generated by means of the Leath algorithm (Aharony and Stauffer, 1992) in order to avoid the presence of finite subclusters and to ensure the connectivity of the solid reactant matrix. The fluid species A diffuses through the product layer with constant diffusivity DAC. As the reaction proceeds, the product layer progressively fills the complex space initially occupied by the solid reactant. Linear intrinsic kinetics, ψ(ξ) ) 1, is assumed. Figure 3 shows the behavior of the total conversion X vs τ, in the diffusion-controlled regime (R ) 0.1, φ2 ) 4), obtained by applying the lattice simulator (continuous lines) in the case of slab pellets characterized by different values of the (statistically) uniform solid reactant fraction . The SC model can still be applied in this case (diffusion-controlled regime) by taking the effective surface Thiele modulus as equal to φ˜ 2 ) φ2/ DΣ. Here DΣ ) DΣ() is the effective diffusivity in the space complementary to the inert matrix (with macroporosity given by the solid reactant fraction ) such that the effective diffusivity of A in C within the product layer, spatially confined by the inert matrix, is given by DACe ) DACDΣ. The values of DΣ() for spatially uncorrelated percolation lattices far from criticality ( g 0.7) have been independently evaluated from the solution of the exittime equation (Giona et al., 1995) and are shown in

Figure 4. Effective diffusivity DΣ vs obtained from the solution of the exit-time equation (dots) (see Giona et al., 1995). (a) Twodimensional spatially uncorrelated pore structures. The continuous line corresponds to the renormalization group result obtained by Shah and Ottino (1986). (b) Two-dimensional Gaussian correlated pore structures. λ(∆η)2 ) 5 × 10-2. (c) Gaussian correlated pore structures. λ(∆η)2 ) 1 × 10-2.

Figure 5. X vs τ as furnished by the lattice simulator (continuous lines, N ) 200) for an initially nonporous slab particle. The fluid species A diffuses through the product layer with a concentrationdependent diffusion coefficient DAC(ξ) ) DACoξ. The intrinsic kinetics is linear, and the phenomenon is diffusion-controlled. Dots represent the solution of the corresponding dynamic SC model. Curve a refers to an initially nonporous pure solid reactant particle, ) 1 (φ˜ 2 ) φ2 ) k(cA0)L/DACo in the corresponding SC model, circles). Curve c refers to an initially nonporous slab particle with solid reactant fraction ) 0.8 (φ˜ 2 ) φ2/DΣ in the corresponding SC model, circles). For both curves, R ) 0.1 and φ2 ) 10. The dotted lines (curves b and d) are the corresponding lattice simulation results for a constant diffusivity DAC ) DACo.

Figure 4 (curve a) together with the renormalization group results obtained by Shah and Ottino (1986). Since the percolation lattice is statistically spatially uncorrelated, the macroporosity of the product layer for each ηc, i.e., for each conversion, is uniform and given by . The accuracy of lattice simulation results is confirmed by the excellent agreement with the solution of the dynamic SC model, eqs 18-20 (dots), with φ˜ 2 ) φ2/DΣ. The accuracy of simulation results has been checked also for low values of R (corresponding qualitatively to gas-solid systems) and for higher values of the Thiele modulus φ. Analogously, simulation results have been quantitalively checked versus the corresponding SC model in the case of the concentration-dependent diffusion coefficient DAC(ξ) ) DACoξ, as shown in Figure 5. In Figure 5, curve a refers to an initially nonporous pure solid reactant particle, ) 1 (φ˜ 2 ) φ2 ) k(cA0)L/DACo in the corresponding SC model, dots). Curve c refers to an initially nonporous slab particle with solid reactant fraction ) 0.8 (φ˜ 2 ) φ2/DΣ in the corresponding SC


in the inert matrix (solid reactant fraction ) 0.8). Curves a and b refer to a Langmuir-Hinshelwood kinetics:

-RA )

kscA 1 , k(cA0) ) ks , ψ(ξ) ) 1 + KcA 1 + Kξ

(24)

for two different values of the K ) KcA0 parameter. Curve c refers to a second-order kinetics:

-RA ) kccA2, ψ(ξ) ) ξ, k(cA0) ) kccA0 Figure 6. X(τ) vs τ furnished by the lattice simulator (continuous lines, N ) 150) for an initially nonporous cylindrical particle formed by the solid reactant B dispersed in an inert matrix for different values of the uniform solid reactant fraction , in the diffusion-controlled regime. The intrinsic kinetics is linear, and the R value (R ) 0.1) corresponds to a liquid-solid system. (a) ) 1.0. (b) ) 0.9. (c) ) 0.8. (d) ) 0.75. dots are the solution of the dynamic SC model (R ) 0.1, φ2 ) 4) with φ˜ 2 ) φ2/DΣ().

(25)

Dots are the solution of the corresponding SC model. The quantitative accuracy of simulation results in the case of nonlinear kinetics is confirmed by the excellent agreement with the SC model in the case of n-order and Langmuir-Hinshelwood kinetics (Figure 7). There are, however, cases for which the results obtained by the simulator deviate from the corresponding SC model, due to the intrinsic inability of the SC model to describe the reaction in the presence of roughening of the reaction interfaces. A typical example is the modified Langmuir-Hinshelwood kinetics, studied by Erk and Dudukovic (1984):

klcA

-RA )

, (1 + KlcA)2

ψ(ξ) )

1 , (1 + Klξ)2 k(cA0) ) kl (26)

Figure 7. X(τ) vs τ furnished by the lattice simulator (continuous lines) in the case of an initially nonporous slab pellet (N ) 200) (solid reactant fraction ) 0.8) for two different nonlinear intrinsic kinetics. Curves a and b refer to a Langmuir-Hinshelwood kinetics with φ2 ) 4.0 and R ) 0.1. (a) K ) 0.1; (b) K ) 1. Curve c refers to a second-order kinetics, φ2 ) 4.0 and R ) 0.1. Dots represent the solution of the corresponding SC model.

model, dots). For both curves, R ) 0.1 and φ2 ) 10. The dotted lines b and d are the corresponding lattice simulation results for a constant diffusivity DAC ) DACo. It is important to stress that the comparison with the continuous SC model does not contain any fitting (adjustable) parameter, since the values of DΣ were evaluated independently and are a characteristic property of the lattice configuration, while the other parameters R, , and φ are fixed and coincide in simulations and continuum solutions. The lattice simulator furnishes accurate results also for cylindrical pellets. Figure 6 shows the same results as Figure 3, but in the case of cylindrical pellets. The satisfactory agreement between simulation results (continuous lines) and the solution of the dynamic SC model (dots), eqs 18-20 with d ) 2 and φ˜ 2 ) φ2/DΣ, confirms the correct choice of local rescaling, eq 19, by the geometric angular factors Zij eq 18, accounting for finite curvature effects. 3.3. Nonlinear Intrinsic Kinetics. The lattice simulator was developed for a generic nonlinear intrinsic kinetics described by the dimensionless functional ψ(ξ). Figure 7 shows the behavior of the total conversion X vs τ in the case of an initially nonporous slab pellet, constituted by the solid reactant B randomly distributed

For this kinetics, it has been shown that, in the diffusion-controlled regime and for negligible external mass-transfer resistance, geometric instabilities of the reaction front occur for slab, cylindrical, and spherical particles, due to the nonmonotonic behavior of the reaction rate -RA as a function of the fluid concentration cA (Erk and Dudukovic, 1984). In particular, nonuniformities in the solid reactant matrix determine nonhomogeneities in the reaction front, and the selfinhibited kinetics causes higher reaction rates at points of deeper penetration. This implies an ever-increasing unevenness of the reaction front, for values of ξ greater than ξm(K ), corresponding to the local maximum of ψ(ξ). In this case the SC approximation no longer holds true. Figure 8 shows the simulation results (continuous lines) for an initially nonporous slab particle (solid reactant fraction ) 0.8) in the presence of a modified Langmuir-Hinshelwood kinetics (Kl ) 4) for several values of the surface Thiele modulus φ. Let ξm(K ) be the value of the dimensionless fluid concentration ξ in correspondence with which the reaction rate -RA attains its maximum value and such that, for ξ > ξm, reaction-front instability may occur and let Xm(φ) be the value of the total conversion in correspondence with which the average fluid concentration on the solid reaction front attains the value ξm and consequently the curve X vs τ exhibits an inflection point. Figure 8 clearly shows that simulation results also exhibit an inflection point which moves toward lower conversion values as the surface Thiele modulus φ increases, as theoretically expected in the diffusioncontrolled regime. The same figure also shows the comparison with the solution of the corresponding dynamic SC model (dotted lines).


ξ ) cA/cA0,

η ) x/L,

θ ) T/T0,

τ ) tkT/L (28)

and the dimensionless parameters

φ2 )

kTL , DAC σ)

Figure 8. X(τ) vs τ for an initially nonporous slab particle (N ) 200) formed by the solid reactant B randomly distributed in an inert matrix (solid reactant fraction ) 0.8) in the presence of a modified Langmuir-Hinshelwood kinetics, eq 26 in the diffusioncontrolled regime (R ) 0.1 and Kl ) 4). (a) φ2 ) 16. (b) φ2 ) 25. (c) φ2 ) 50. (d) ) 0.8. Dotted curves are the solution of the corresponding dynamic SC model, eqs 21-23 and eq 26 with φ˜ 2 ) φ2/DΣ.

It should be observed that, for each value of φ, the SC model deviates significantly from simulation results for X < Xm and, after the inflection point is reached, the SC model and the simulation results proceed parallel to one another. Indeed, for X < Xm (i.e., for ξ > ξm) the disordered nature of the solid matrix induces reaction-front instability, and the SC model does not correctly reproduce the macroscopic properties of the system. However, for X > Xm the modified Langmuir-Hinshelwood kinetics behaves as a monotonically increasing function of fluid concentration, and reaction-front instabilities are progressively reabsorbed. Under these conditions, SC approximation is again valid. As confirmation of this analysis, Figure 8 clearly shows that the divergence between the SC model and simulation results decreases as the Thiele modulus φ increases, i.e., as the characteristic length scale of the instability zone (occurring for X < Xm) decreases. The detailed scaling analysis of the influence of nonmonotonicity in the presence of disordered solid mixtures (inert + reactant) lies beyond the scope of this article and will be discussed elsewhere. It is, however, important to observe that the case of the modified Langmuir-Hinshelwood kinetics is a typical example in which the lattice simulator shows all its usefulness, since it enables us to study the influence of disorder on the evolution of reaction both globally (conversion/time plots) and microscopically (roughnening of the reaction front). Similar data cannot be obtained by means of continuum models. 3.4. Nonisothermal Intrinsic Kinetics. The way the diffusion step is implemented and the validity of the reaction step for generic nonlinear kinetics lead us to address nonisothermal intrinsic kinetics of the form

-RA ) kTcA exp(-E/RT)

(27)

The case of a nonlinear dependence of -RA on cA can be worked out, as discussed in Section 3.3. For the sake of simplicity, let us consider the case of a pure solid reactant particle. Let T, (Fcp)P, and KP be respectively the temperature field, the volumetric heat capacity, and thermal diffusivity in the product layer. Let us further introduce the following dimensionless variables:

φT2 )

E , RT

kTL , KP

γ)

R)

bcA0 FB

cA0(-∆H) (Fcp)PT0

(29)

We shall consider the simplified case of a small (negligible) volumetric heat capacity of the solid reactant (Fcp)B with respect to (Fcp)P, i.e., ζB ) (Fcp)B/(Fcp)P f 0. While this assumption simplifies the formal treatment of the problem, the general expressions for the discrete energy balance are also given. In the case of nonisothermal kinetics, eq 9 attains the form

∆η ) Rξ exp(-σ/θ) ∆τ

(30)

and the dimensionless time interval ∆τ′ can be defined as ∆τ′ ) ∆η/R, corresponding to the dimensionless time scale required for the advancement of the local reaction front by the lattice length scale ∆η at the maximum fluid concentration ξ ) 1 and at infinite temperature θ ) ∞. In the diffusive step, both the concentration and the temperature fields ξij and θij within the product layer ((i, j) ∈ P) are updated according to the diffusive dynamics. Indeed, the assumption that ζB f 0 implies that the temperature increment of the unreacted solid may be overlooked in heat transfer. For a constant diffusivity DAC, the diffusive dynamics of the fluid species A remains unchanged with respect to the isothermal case (eq 12). For a temperaturedependent diffusion coefficient, DAC(θ) ) DACTg(θ), the algorithm introduced for a concentration-dependent diffusion coefficient (eq 13) can be extended in a straighforward way by simply replacing f(ξ) and φ2 ) k(cA0)L/DAC0 respectively with g(θ) and φ2 ) kTL/DACT. In the same way, the temperature field is updated in the diffusive step according to the relation

θij(n+1) ) θij(n) +

ωT 4

∑

χhk[θhk(n) - θij(n)],

(h,k)∈I(i,j)

n ) 0, ..., NT (31) to be repeated NT times, where NT ) [4(∆τs/∆η2)/φT2] and ωT ) 4(∆τs/∆η2)/(NTφT2). Of course θij(0) ) θij(τ) and θij(NT) ) θij(τ + ∆τs/2). Equation 31 holds only in the case of ζB f 0. In the more general case (ζB > 0), the temperature field should be defined for each lattice site (i, j), and the heat conduction algorithm attains the following form:

θij(n+1) ) θij(n) + ωT

∑

2

KijKhkζhk

[θhk(n) - θij(n)] (32) K K ζ + K ζ hk hk 4K ζ(h,k)∈I(i,j) P ij ij

where Kij ) KP if (i, j) ∈ P, Kij ) KB if (i, j) ∈ B, ζij ) ζB if (i, j) ∈ B, and ζij ) 1 if (i, j) ∈ P. For each lattice site, the recursive relation eq 32 is repeated NT times, where


ξ ) 1 at η ) 1,

θ ) 1 at η ) 1,

∂ξ ) (φ2)ξ exp(-σ/θ) ∂η at η ) ηc (38) ∂θ ) -(φT2)γξ exp(-σ/θ) ∂η at η ) ηc (39)

dηc ) -R(ξ exp(-σ/θ))|η)ηc dτ

Figure 9. X(τ) vs τ furnished by the lattice simulator (continuous lines, N ) 200) for an initially nonporous slab particle formed by the pure solid reactant in the case of a nonisothermal intrinsic kinetics (eq 27). (a) σ ) 4.5, ψ ) 5.0, φ2 ) 0.01, φT2 ) 4, and R ) 4.0. (b) σ ) 4.5, ψ ) 5.0, φ2 ) 1, φT2 ) 4, R ) 1.0. (c) σ ) 3.0, ψ ) 4, φ2 ) 10, φT2 ) 3, and R ) 0.5. Continuous lines are the solution of the corresponding dynamic SC model eqs 36-40.

NT ) [(∆τsKζ/∆η2)/φT2], K ) max{KP, 1} and ζ ) max{1, ζB}. Equation 32 reduces to eq 31 for ζB f 0. The origin of eq 32 is discussed in the Appendix by deriving the discrete version of the thermal balance equation, in the case of position-dependent volumetric heat capacity and thermal conductivity. In the reacting step, the solid reactant fraction qij in each site (i, j) ∈ B and the concentration field ξhk for (h, k) ∈ I(i,j) ∈ P are updated according to eqs 15 and 16 with

∆qij ) 〈ξij〉 exp(-σ/〈θij〉)

(33)

where the symbol 〈 〉 stands for the local average 〈θij〉 ) (1/νij)Σ(h,k)∈I(i,j)χhkθhk, with νij ) Σ(h,k)∈I(i,j)χhk. Due to the nonisothermal nature of the reaction, heat generation occurs at the fluid-solid interface. Consequently, the temperature field should be updated within the reacting step to account for it. In the case ζB f 0, the heat generated is redistributed only between those sites (h, k) adjacent to the solid reactant site (i, j) belonging to the product layer, so that

(

θhk f θhk 1 +

f ∆qijγ R νij〈θij〉

)

(h, k) ∈ P

(34)

Equation 34 is analogous to eq 16 for updating the fluid concentration field due to the reaction. In the case ζB > 0, the enthalpic balance leads to the following updated rule:

(

ζhkθhk f ζhkθhk 1 +

∆qijγ f R (p + 1)〈ζijθij〉

)

(35)

where p ) 4 is the number of nearest-neighboring sites for a regular two-dimensional square lattice. Figure 9 shows the comparison between simulation results (dots) and the solution of the corresponding dynamic nonisothermal SC model (continuous line) for ζB f 0:

1 ∂2ξ ∂ξ ) 2 2, ∂τ φ ∂η

ηc e η e 1

(36)

∂θ 1 ∂2θ ) 2 2, ∂τ φ ∂η

ηc e η e 1

(37)

T

(40)

The agreement is fully satisfactory. The different curves refer to three different operating conditions, two of which produce a sudden ignition phenomenon, as can be observed both in lattice simulations and in the continuum solution. 4. Further Applications The examples discussed in Section 3 refer to initially nonporous particles formed by one solid reactant dispersed in an inert matrix in the presence of only one fluid reactant. The algorithm can easily be extended to more complex cases, such as (1) multicomponent diffusion, (2) n solid reactant species, and (3) initially porous particles. This section briefly analyzes these cases. 4.1. Multicomponent Diffusion. Let us consider the following fluid-solid noncatalytic reaction:

a1A1(f) + a2A2(f) + B(s) f cC(f) + dD(s)

(41)

where, to simplify the notation, we assume c ) a1 + a2, i.e., equimolecular diffusion of the fluid species in the product layer. The initially nonporous pellet is formed by the solid reactant B dispersed in an inert matrix and the intrinsic kinetics is given by

( )

-RB ) km(cA10,cA20)cA1cA2ψ

cA1 cA2 , cA10 cA20

(42)

where cA1 and cA2 are the concentrations of the fluid species in the product layer, and cA10 and cA20 are the external fluid concentrations. Let us introduce the dimensionless variables

x η) , L

τ)

tkmcA10 , L

ξ1 )

cA1 0

c A1

,

ξ2 )

cA2 cA20

(43)

and the dimensionless parameters

kmcA10L , φ ) D 2

cA20 R) , FB

β)

cA20 cA10

(44)

where D ) max{D1m, D2m} and the diffusivities D1m and D2m of both A1 and A2 in the fluid mixture are assumed to be constant. In this case, the dimensionless time interval ∆τ′ is given by ∆τ′ ) ∆η/(Rψ(1,1)). At each sampling instant, ∆τs ) f∆τ′, the concentration fields ξij1 and ξij2 for (i, j) ∈ P and qij for (i, j) ∈ B are updated according to the previously defined timesplitting algorithm. In the diffusion step, the diffusive dynamics of the fluid species in the product layer is given by the


following algorithm:

ξijl(n+1) ) ξijl(n) + ωDDlm

∑

4D

χhk[ξhkl(n) - ξijl(n)]

l ) 1, 2 (45)

(h,k)∈I(i,j)

to be repeated Nd ) [4∆τs/(∆η2φ2)] times, with ωD ) 4∆τs/ (Nd∆η2φ2). In the reaction step, the solid fraction ∆qij at site (i, j) converted by the reaction is given by

〈ξij 〉〈ξij 〉ψ(〈ξij 〉,〈ξij 〉) 1

∆qij )

2

1

(46)

) )

f a1β∆qij ξhk1 f ξhk1(τ) 1 - χhk R 〈ξ 1〉ν ij ij

(47)

f a2∆qij ξhk2 f ξhk2(τ) 1 - χhk R 〈ξ 2〉ν

(48)

ij

ij

The case of concentration-dependent Fickian diffusion coefficients, which arise if one adopts the general (nonlinearized) expression for the fluxes (StefanMaxwell equations) or are induced by thermodynamic nonidealities in the fluid phase (Taylor and Krishna, 1993), can be treated by merging eq 45 and eq 13. 4.2. Solid Reactant Mixtures. Let us consider the case of a solid reactant mixture characterized by two species, B1 and B2, undergoing the parallel/competitive reactions

A(f) + b1B1(s) f C1(f) + dD(s) A(f) + b2B2(s) f C2(f) + dD(s)

(49)

The fluid species A diffuses through the product layer and reacts on the solid boundary, with the two intrinsic kinetics

-RA1 ) k1(cA0)ψ1(cA/cA0) -RA2 ) k2(cA0)ψ2(cA/cA0)

(50)

By defining the dimensionless time τ with respect to the maximum reaction rate k ) max{k1(cA0), k2(cA0)}, τ ) tk/L, the characteristic time interval ∆τ′ is given by

{

}

∆η ∆η ∆η ) , 2 2 1 R ψ (1) R ψ (1) Rψ(1) 1

Rl〈ξij〉ψl(〈ξij〉)

(52)

Rψ(1)

In order to consume f∆qijFB∆x2 mol of the solid reactant site (i, j), by enforcing the stoichiometry, f∆qijFB∆x2a1 moles of the fluid species A1 and f ∆qijFB∆x2ba2 mol of the fluid species A2 in the neighborhood of (i, j) are needed. In dimensionless terms, this means that the average fluid concentrations 〈ξij1〉 and 〈ξij2〉 in the neighborhood of the reactant site must be updated as follows 〈ξij1〉 f 〈ξij1〉 - f∆qija1β/R and 〈ξij2〉 f 〈ξij2〉 f∆qija2/R. This is obtained by updating the concentration fields of the fluid species A1 and A2 at fluid sites (h, k) ∈ I(i,j) belonging to the product layer as follows:

∆τ′ ) min

∆qijl )

1

ψ(1,1)

( (

The diffusion step is analogous to the case of a single solid reactant. In the presence of two solid reactant species B1 and B2, we define two continuous concentration fields qij1 and qhk2 representing respectively the unreacted fractions of solid reactants B1 and B2 at sites (i, j) ∈ B1 and (h, k) ∈ B2, where Bi, (i ) 1, 2) is the set of lattice sites formed by the ith solid reactant. In the reaction step, the concentration fields {qij1}, {qij2}, and {ξij} are updated as

(51)

where Rl ) blcA0/Fl and Rψ(1) ) max{Rlψl(1)}, for l ) 1, 2.

qijl(τ + ∆τs) ) qijl(τ) - f∆qijl

(

(i, j) ∈ Bl,

)

l f ∆qij ξhk f ξhk(τ) 1 - χhk l R 〈ξij〉νij

l ) 1, 2 (53)

(h, k) ∈ I(i,j) (54)

The value of f should be chosen in such a way that max{f/R1, f/R2} < 1. 4.3. Initially Porous Particles. In many fluidsolid systems, the solid matrix is initially porous, allowing diffusion and a reaction to occur simultaneously throughout the entire pore space, which evolves in time. Therefore, the reaction occurs in a broad zone rather than at a sharp boundary, and SC approximation is no longer applicable. In many instances, a chemical reaction consumes the porous medium, leading eventually to the total consumption and disappearence of the solid matrix. In other cases, reaction products may be deposited on the solid matrix, resulting in decreased reactivity and pore blockage (Sotirchos and Yu, 1985; Sotirchos and Zarkanitis, 1993; Yu and Sotirchos, 1987; Ramachadran and Smith, 1977; Salles et al., 1993; Bekri et al., 1995). Both these phenomenona can be addressed by means of the lattice simulator discussed in the previous sections for initially nonporous particles with no significant modifications since, from the algorithmic point of view, it makes no difference whether the pore space is initially present or progressively created by the consumption of the solid reactant. As an example, Figure 10 shows the dissolution process occurring on a two-dimensional initially porous particle formed by a pure solid reactant B (in black), in the case where the initial pore space is represented by a diffusion limited aggregate (DLA) (Havlin and BenAvraham, 1987). The intrinsic kinetics is linear in the fluid concentration, and dissolution occurs in the reaction-controlled regime (φ2 ) 10-3), i.e., the concentration of fluid reactant is essentially uniform throughout the entire pore space and equal to the maximum external fluid concentration. For this system, the fragmentation process is very important due to the finitely ramified and loopless nature of the DLA cluster and to the deep penetration of the diffusion front (Sahimi and Tsotsis, 1987, 1988; Kerstein and Edwards, 1987; Cai et al., 1991; Gyure and Edwards, 1992). Two competing phenomena determine the evolution of the pore network: (1) pore growth and (2) coalescence of neighboring pores. Pore enlargement is the dominant mechanism at the beginning of the process. As the pores get larger, they start to coalesce with adjacent pores, thus causing the surface area to decrease. Figure 11A shows the number Nc of finite clusters (fragments) in which the solid structure disintegrates


Figure 10. Dissolution process occurring in a two-dimensional initially porous particle of pure solid reactant B (in black, N ) 200) in the case where the initial pore space is represented by a DLA structure. The intrinsic kinetics is linear in the fluid concentration, and dissolution occurs in the reaction-controlled regime (φ2 ) 10-3). (A) Initial configuration of the pore space (X ) 0); (B) X ) 0.25; (C) X ) 0.5.

Figure 11. (A) Number Nc of finite clusters (fragments) in which the solid structure of Figure 10 disintegrates as a function of the total conversion X. (B) Behavior of the total conversion X as a function of time τ for the dissolution process depicted in Figure 10.

as a function of the total conversion and Figure 11B the behavior of the total conversion X as a function of time τ. At low conversion, the number of fragments increases and reaches a maximum. In this case reaction occurs at the smallest length scales and connects together many dangling ends of the initially loopless pore structure, thus generating an increasing number of small clusters. Correspondingly, the total conversion X increases in linear fashion with time up to a given point Xm, at which the curve Nc(X) reaches its local maximum. For X > Xm, the global conversion curve starts to inflect and the total number of clusters decreases. This derives from the fact that maximum fragmentation has been reached and the reaction shrinks each fragment progressively. Dissolution simulations in the diffusioncontrolled regime can be performed in exactly the same way. 5. Influence of Spatial Correlations Random percolation networks are an important tool for modeling transport processes in disordered systems such as porous media, polymers, and gels. In the study of percolation models, spatial disorder is usually assumed to be uncorrelated, i.e., the probability of any site being occupied is independent of the occupancy of other sites. However, natural systems also exhibit some degree of spatial correlation, and correlated percolation models prove useful in application to porous media and transport because they mimic the structure of real materials better than the customary uncorrelated percolation schemes. The generation of correlated percolation lattices can be regarded as an application of the reconstruction of

porous media. This topic has been extensively studied by Joshi (1974), Quiblier (1984), and Adler et al. (1990). An analytic solution of the inverse reconstruction problem has recently been proposed by Giona and Adrover (1996) for isotropic and homogeneous porous media and further extended to the generation of nonhomogeneous porous matrices with a prescribed pore-pore correlation function and position-dependent porosity (Adrover and Giona, 1997). In this section we analyze the influence of the spatial correlation of solid reactant distribution on fluid-solid noncatalytic reactions, focusing on initially nonporous particles with an assigned solid reactant distribution enclosed within an inert matrix (Sohn and Xia, 1986, 1987). Let us consider a thin section of an initially nonporous pellet consisting of solid reactant B and inert solid I, described by means of the solid reactant characteristic function, χB(r) ) 1 if r ∈ B and χB ) 0 otherwise, r being the lattice-position vector. The spatial correlation properties of the solid reactant matrix are described by means of the second-order normalized correlation function

C2χB(x) )

〈(χB(r) - )(χB(r + x) - )〉 - 2

(55)

where ) 〈χB(r)〉 is the uniform solid reactant fraction and 〈 〉 stands for spatial average. Under the assumption of isotropy, C2χB(x) is solely a function of x ) |x| and a generic cross-section of the material is representative of the entire structure (Adler et al., 1990). Figure 12 shows two cylindrical pellets (solid reactant in black, inert matrix in white) characterized by the same uniform solid reactant fraction, ) 0.8, and by a Gaussian decay in the spatial correlation function C2χB ) exp(-λx2), for two different values of the correlation length Lc ) 2/(λ)1/2. For isotropic, nonhomogeneous porous pellets characterized by a radially symmetric solid reactant distribution (r), the second-order correlation function C2χB(x) attains the form

C2χB(x) )

〈(χB(r) - (r))(χB(r + x) - (|r + x|))〉 〈(r)〉 - 〈2(r)〉

(56)

By controlling the spatial correlation properties and the radial distribution of the solid reactant, it is possible


Figure 12. Representation of cylindrical pellets (solid reactant in white, inert matrix in black) characterized by the same uniform solid reactant fraction ) 0.8 and Gaussian decay of the spatial correlation function C2χB(x) ) exp(-λx2) for two different values of the correlation length Lc ) 2/(λ)1/2. (A) λ(∆η)2 ) 0.5; (B) λ(∆η)2 ) 0.05.

Figure 13. X(τ) vs τ furnished by lattice simulations, in the case of a linear intrinsic kinetics (eqs 1-2), in cylindrical pellets (N ) 150) characterized by a uniform solid reactant fraction and by a Gaussian decay of the spatial correlation function (continuous lines, λ(∆η)2 ) 0.05). The simulation parameters are φ2 ) 4 and R ) 0.1. (a) ) 0.8. (b) ) 0.7. Dotted lines represent lattice simulation results for the corresponding uncorrelated structures, with the same . Dots represent the solution of the corresponding SC model (eqs 21-23) with φ˜ 2 ) φ2/DΣ.

Figure 14. X(τ) vs τ furnished by lattice simulations in the case of a linear intrinsic kinetics (eqs 1-2) in cylindrical pellets (N ) 150) characterized by a linear solid reactant distribution (η) ) 0 + (R - 0)η and by a Gaussian decay of the spatial correlation function (continuous lines, λ(∆η)2 ) 0.05). The simulation parameters are φ2 ) 1 and R ) 0.1. (a) 0 ) 0.7 and R ) 1.0. (b) 0 ) 1.0 and R ) 0.7. Dotted lines represents lattice simulation results for the corresponding uncorrelated structures, with the same (η).

to generate solid pellets with complex geometry and connectivity. Indeed, the pore space (or product layer space) which originates from the solid reactant consumption can exhibit a distributed porosity and a wide pore-size distribution. Let us first consider the case of uniform solid reactant distribution, i.e., (r) ) ) constant. Figure 13 shows the behavior of X vs τ obtained by lattice simulations in the case of a linear intrinsic kinetics (eqs 1 and 2) in cylindrical pellets characterized by a uniform solid reactant fraction and by a Gaussian decay of the spatial correlation function (continuous lines). Dotted lines represent lattice simulation results for the corresponding uncorrelated structures, with the same , and dots are the solution of the corresponding SC model (eqs 21-23) with φ˜ 2 ) φ2/DΣ. The behavior of DΣ vs for different values of λ is shown in Figure 4 (these data were obtained by solving the exit-time equation; see Giona et al. (1996)). For a given , the effective diffusivity increases with the increasing of the correlation length, but the effects of correlation on DΣ are significant only for low values of (i.e., approaching the critical threshold c = 0.593 for d ) 2). The good level of agreement between lattice simulation results and the SC model (see Figure 13) in the

diffusion-controlled regime (φ2 ) 4.0, R ) 0.1) once again confirms the quantitative accuracy of the lattice simulations and the possibility of describing the effects of the complex and disordered nature of the solid reactant matrix by simply introducing the effective diffusivity of the fluid species in the space surrounded by the inert matrix and progressively filled by the product layer. Many authors have analyzed the influence of the solid reactant distribution on fluid-solid noncatalytic reactions overlooking the influence of spatial correlation properties (see Sohn and Xia (1986, 1987) and references therein). Here we study the influence of correlations on fluidsolid noncatalytic reactions for initially nonporous particles characterized by a linear solid reactant distribution (Dudukovic, 1984) (η) ) 0 + (R - 0)η and by a Gaussian decay of the correlation function. Figure 14 shows the behavior of X(τ) vs τ for two correlated initially nonporous particles characterized by the same value of λ and by two different linear solid distribution functions, compared with the corresponding conversion/time curves of uncorrelated structures (dotted curves). As might be expected (since effective diffusivity increases with correlation length), the total conversion time is smaller for correlated structures than


Figure 15. φ˜ L2(η)/φ2 vs X for Gaussian correlated (solid lines, λ(∆η)2 ) 0.05) and uncorrelated structures (dotted lines) with linear solid fraction distribution. (a) 0 ) 0.7 and R ) 1.0. (b) 0 ) 1.0 and R ) 0.7.

for uncorrelated ones. Curve a refers to a particle with an increasing solid fraction with the dimensionless radius η (0 ) 0.7, R ) 1.0). In this case, the effects of correlations are weak and appear only at high conversions. Curve b refers to a particle possessing a solid fraction decreasing with η, (0 ) 1.0, R ) 0.7). Spatial correlation effects are significant also at low conversion. It is evident that the influence of correlations is significant only for those structures which possess lower solid fractions in contact with higher values of fluid concentration. This result can be easily explained by analyzing the behavior of the local effective Thiele modulus, defined as 2

φ˜ L (ηc) )

φ2(1 - ηc)(ηc) DΣ((ηc))

(57)

as a function of the conversion X(ηc). Figure 15 shows the the behavior of the effective Thiele modulus, defined by eq 57, of correlated structures compared with that of uncorrelated ones. In the case of a solid fraction increasing with η, the curves φ˜ L2(η)/φ2 vs X(η) for correlated and uncorrelated structures practically coincide up to high values of conversion (Figure 15, curve a). In the case of solid fractions decreasing with η, the curves φ˜ L2(η)/φ2 vs X(η) associated with correlated and uncorrelated structures diverge significantly, especially for intermediate values of the conversion 0.1 e X e 0.6. This effect (deviation) shows up in the X-τ curves, which initially deviate appreciably and then proceed practically parallel to each other. To sum up, the influence of porosity and spatial correlations on fluid-solid noncatalytic reactions involving initially nonporous particles (and monotonic kinetics) can be completely described by means of the effective medium diffusivity depending on both and λ. The effect of spatial correlations on the global kinetics is significant only at low solid fractions and for particles presenting low solid fractions in contact with higher values of the fluid concentration. 6. Concluding Remarks The lattice simulator developed in this article proves to be applicable for arbitrary nonlinear kinetics also in nonisothermal conditions, as shown by the excellent agreement between simulation results and the SC model analyzed in Section 3. Indeed, none of the comparisons between simulation and continuum models contain

adjustable parameters since the effective transport parameters entering into the SC model are obtained independently of the reaction considered and constitute intrinsic properties of the lattice structure. From the thorough analysis developed in Section 3, it can be observed that the lattice simulator proposed in this article overcomes all the limitations of randomwalk and “stick-and-react” models. The definition of a continuum field at each solid reactant site enables us to control and reduce the discretization approximations, which can be crude in all the lattice simulations based on the removal of the entire lattice site once the reactive event has occurred. The diffusive step can easily account for nonlinearities in the Fickian diffusion coefficients (Section 3.2), frequently occurring as a consequence of multicomponent transport and thermodynamic nonidealities (Taylor and Krishna, 1993). The lattice simulator proposed is strictly deterministic. It can be viewed as a mean-field method to approach fluid-solid reaction in complex and disordered matrices, and no fluctuations are present in either the diffusive step or the reactive dynamics. [It should be noted that randomness (thermal fluctuations) can easily be included within the reactive step by making the update of the concentration field a stochastic event while still enforcing the stoichiometry of the reaction.] To some extent, the lattice simulator can be regarded as a computational laboratory by means of which it is possible to study the relationships existing between structural properties of granular pellets, their evolution in time due to reaction, and the global kinetics. For example, in the case of initially nonporous particles, the simulator can be applied to achieve simultaneous analysis of the overall kinetics and the form and shape (roughness, groove formations, etc.) of the evolving reaction interface. In this way, the effects of the structural heterogeneity (randomness in the distribution of the solid reactant within an inert matrix) on the morphological properties of the evolving front can be studied in detail. Similarly, in the case of porous solids, the lattice simulator proposed makes it possible to analyze at a microscopic and statistical level the structural features of the evolving structure, such as the accessible surface area, the pore-pore correlation function, etc. Such data are invaluable in the formulation of accurate continuum models and in the definition of more detailed mechanistic descriptions of the influence of structural and textural features characterizing porous structures on overall reaction kinetics. Let us briefly address the computational issues. The computer time and memory required to implement the simulator depend strongly on the operating conditions since the integers Nd and NT entering into the diffusive step depend on φ2 and the parameter f depends on R (due to stability constraints). In all the simulations performed throughout this article, lattice size is taken as equal to 400 × 400, i.e., N ) 200, and the computer time required to achieve complete conversion (on a RISC workstation) was under 10 h in the worst case. It should be observed that the structure of the diffusive step and the reaction step is particularly suitable for parallel implementation, and the computer time can therefore be sharply reduced. This article mainly analyzes moderately disordered structures. The discussion developed in the Appendix also shows that the algorithm can be applied in the presence of highly singular structures such as fractals. The article analyzes in some detail the effect of spatial


correlation properties on the distribution of the solid reactant. The transport (and consequently reaction) properties of short-range correlated structures differ in a significant way from those of pure random distributions. The results obtained indicate that continuum modeling should also take into account the effects of correlations. From the experimental point of view, correlation properties can be measured directly by analyzing digitized images of two-dimensional thin sections or indirectly by means of scattering experiments (since the structure factor gives the Fourier transform of the pair correlation function). Levitz (1996) has recently proposed some interesting approaches to relate the correlation properties of porous structure and scattering data. Our analysis has focused almost exclusively on initially nonporous particles. For monotonic reaction, the influence of disorder and correlation shows itself exclusively through a modified effective diffusivity. The situation is slightly different in the case of nonmonotonic reactions (such as the modified Langmuir-Hinshelwood kinetics). In this case, small heterogeneities and the local clustering of the reactant may cause the formation of grooves in the reaction interface which propagate faster and thus induce a progessive roughening of the interface for fluid concentrations less than the critical concentration corresponding to the local maximum of the intrinsic rate function. In accordance with the intrinsically methodological purpose of this article, we have shown that generic fluid-solid reactions can be studied with accuracy on arbitrary complex lattice configurations. The method proposed can be viewed as an algorithm to solve movingboundary problems in complex geometries and opens up interesting perspectives in the study of all the phenomenologies associated with reactions on initially porous particles, in the presence of reliable models for multicomponent transport and reactions starting from a strict geometric lattice representation of the solid particles. Acknowledgment Authors thank Alessandro Galassini for useful discussion. This work was partially supported by MURST 40% grants. Notation B:

set of sites belonging to the solid reactant matrix

C:

set of sites belonging to the inert matrix

C2χB:

normalized pore-pore correlation function of the d-dimensional reconstructed porous structure

d:

Euclidean-space dimension

df:

fractal dimension

d s:

spectral dimension

d w:

walk dimension

DAC: e:

diffusivity of A in C in the product layer

DAC

DACDΣ effective diffusivity of A in C in the product layer, confined in the space complementary to the inert matrix

K:

thermal conductivity

Mt/M∞:

fractional mass uptake

N:

lattice size [lattice units]

R2(t):

mean square displacement

P:

set of sites belonging to the product layer

qij:

unreacted fraction of solid reactant at site (i, j)

Zij:

cylindrical geometrical factor

Greek Letters R:

bcA0/FB

γ:

(cA0(-∆H))/((Fcp)PT0), eq 29

∆τ′:

dimensionless time interval

∆τs:

dimensionless sampling time interval

:

solid reactant fraction in the pellet

ζB:

(Fcp)B/(Fcp)P

η:

dimensionless spatial coordinate

ηc:

dimensionless reaction front position

θ:

dimensionless temperature

λ:

Gaussian decay parameter

ξ:

dimensionless fluid concentration

FB:

molar density of pure solid reactant

Fcp:

volumetric heat capacity

σ:

E/(RT), eq 29

τ:

dimensionless time

φ:

surface Thiele modulus, φ2 ) k(cA0)L/DAC

φT:

kTL/KP, eq 29

φ˜ L:

local surface Thiele modulus, eq 57

˜2: φ

φ2/DΣ

χij:

characteristic function

ψ:

dimensionless functional accounting for nonlinearities in the intrisinc kinetics

Superscript and Subscript (n):

order of iteration

ij:

lattice-site coordinates

Special Symbol [x]

smallest integer exceeding x

Appendix This appendix briefly analyzes the numerical method applied in the diffusion step on disordered matrices. The method can be regarded either as a finite-difference approach or as a finite-element approach, as will be explained below. Let us first consider a porous structure (pore networks + solid matrix) and a diffusional dynamics on it, i.e.,

∂c(x,t) ) ∇(D∇c(x,t)) ∂t

(58)

equipped with the zero-flux boundary condition at the solid boundary. D can be a constant, or a positiondependent D(x) or a concentration-dependent diffusion coefficient D(c(x,t))


Figure 17. (A) R2(τ) vs τ ) tD/L2 for the Vicsek fractal (order of iteration n ) 6, N ) 729). The straight line in the log-log plot admits a slope of 2/dw with dw ) df + 1 ) log 5/log 3 + 1. (B) M(τ)/M∞ for the Vicsek fractal fed by the external boundary sites (order of iteration n ) 6). The straight line in the log-log plot admits a slope of df/dw ) df/(df + 1) ) 0.59432.

Figure 16. Vicsek fractal.

The lattice structure is fully described by the characteristic function of its pore network Π, χijΠ:

χijΠ )

{

1 0

if (i, j) ∈ Π otherwise

(59)

In the case of a constant diffusion coefficient, the lattice representation of eq 58 attains the form

cij(t + ∆t) ) cij(t) + ∆t∇d2c|(i,j) 2

∇d c|(i,j) )

1 ∆x2

(60)

∑

Π

χ(h,k) (chk - cij)

(h,k)∈I(i,j)

where cij is the concentration of the transported entity at site (i, j), ∆x the lattice spacing, ∆t the discrete time interval, and ∇d2c|(i,j) the discrete Laplacian operator at site (i, j) ∈ Π. Since the flux from site (h, k) to site (i, j) ∈ Π, J(h,k)f(i,j) is given by

J(h,k)f(i,j) )

D χ Π(c - cij) ∆x hk hk

(61)

it is easy to see that if (h, k) does not belong to Π, zeroflux boundary conditions are enforced automatically. Equation 60 can be viewed as a finite-difference representation of the Laplacian over the lattice quadrature points, or as a finite-element representation in which sites (i, j) ∈ Π are the elements considered, and transport occurs only if two sites (i, j) and (h, k) are connected by a bond (i.e., belong to Π and are nearest neighbors of each other in the prescribed lattice topology). It is important to stress that this approach is the numerical counterpart of the discretization of the balance equations applied in the exact renormalization of transport equations on complex graphs, such as finitely ramified fractals (Giona et al., 1996a-c). The worst numerical case in the application of the algorithm is given by highly singular (and in the limit nondifferentiable) pore-network structures. To give an example of these extreme conditions, we have chosen a finitely ramified fractal, the Vicsek fractal, depicted in Figure 16. Its fractal dimension is given by df ) log 5/log 3 ) 1.465, and since it is a loopless structure its walk dimension dw is df + 1 (Havlin and Ben-Avraham, 1987). Consequently, its spectral dimension is given by ds ) 2df/dw ) log 25/log 15 ) 1.189 (Havlin and BenAvraham, 1987). Two types of boundary-value problems

are considered. The first case is the free evolution of the diffusion front starting from a Dirac’s δ distribution initially located in the center of symmetry of the structure r ) 0. No boundary conditions are assumed, and the simulation stops once the diffusion front reaches the external sites. The characteristic quantity is given by the mean square displacement R2(t) ) ∫|r|2c(r,t) dr, which scales as a function of time as R2(t) ∼ t2/dw. The second case considered is a sorption experiment, starting from an initially empty configuration and with feed from the (four) external boundary sites. In this experiment, the basic quantity is the fractional mass uptake M(t)/M∞, representing the fraction of solute entering the structure up to time t with respect to the saturation value. Exact renormalization analysis (Giona et al., 1996b) shows that in this case M(t)/M∞ ∼ tds/2. Figure 17A and 17B compares the simulation data and the theoretical predictions. Simulations were performed on a Vicsek lattice at iteration n ) 6, corresponding to a 729 × 729 two-dimensional lattice structure. As can be observed, the agreement between simulation and theory is fully satisfactory. This example shows the validity of the simple numerical algorithm applied in the diffusion step also for highly singular porous structures. Similar algorithms have been derived in the case of position-dependent or concentration-dependent diffusion coefficients. Let us first consider the case D(x) ) Dof(x). In this case, the lattice representation of eq 58 attains the form

cij(t + ∆t) ) cij(t) + Do

(

)

1 1 1 1 χ(h,k)ΠFhk,ij(chk - cij) ) + Fhk,ij 2 fhk fij ∆x2(h,k)∈I(i,j) (62) ∆t

∑

where Dofij is the diffusion coefficient at site (i, j) and DoFhk,ij is the average diffusivity over the two adjacent sites (i, j) and (h, k). In a similar way, in the case of a concentrationdependent diffusion coefficient D(x) ) Dof(c(x,t)) eq 58 attains the form

cij(t + ∆t) ) cij(t) + 2f(chk(t)) f(cij(t)) ∆t (chk - cij) (63) χ(h,k)Π Do f(chk(t)) + f(cij(t)) ∆x2(h,k)∈I(i,j)

∑


The algorithm, eq 63, is slightly different if one considers the thermal diffusion equation, in the case where the volumetric heat capacity Fcp and the thermal conductivity k are both functions of the position, i.e.,

Fcp(x)

∂T(x, t) ) ∇(k(x)∇T(x,t)) ∂t

(64)

The discretized version of eq 64 attains the form

Tij(t + ∆t) ) Tij(t) + ∆t 1

∑

2

∆x (Fcp)ij (h,k)∈I(i,j) 1 T

Fhk,ij

)

(

χ(h,k)ΠFThk,ij(Thk - Tij)

)

1 1 1 + 2 khk kij

(65)

which, by introducing the thermal diffusivities Kij ) kij/ (Fcp)ij and Khk ) khk/(Fcp)hk, attains the form

Tij(t + ∆t) ) Tij(t) + 2KijKhk(Fcp)hk χ(h,k)Π (Thk - Tij) (66) Kij(Fcp)ij + Kjk(Fcp)hk ∆x2(h,k)∈I(i,j) ∆t

∑

Literature Cited Adler, P. M.; Jacquin, C. G.; Quiblier, J. A. Flow in simulated porous media. Int. J. Multiphase Flow 1990, 16, 691. Adrover, A.; Giona, M. A Predictive Model for Permeability of Correlated Porous Media. Chem. Eng. J. 1996, 64, 7. Adrover, A.; Giona, M. Reconstruction of Non-Homogeneous Porous Media. Ind. Eng. Chem. Res. 1997, 36, 5010. Aharony, A.; Stauffer, D. Introduction to Percolation Theory; Taylor & Francis: London, 1992. Bekri, S.; Thovert, J. F.; Adler, P. M. Dissolution of porous media. Chem. Eng. Sci. 1995, 50, 2765. Cai, M.; Edwards, B. F.; Han, H. Exact and asymptotic scaling solutions for fragmentation with mass loss. Phys. Rev. A 1991, 43, 656. Crank J. The Mathematics of Diffusion; Claredon Press: Oxford, 1967. Dudukovic, M. P. Reaction of Particles with Nonuniform Distribution of Solid Reactant. The Shrinking Core Model. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 330. Erk, H. F.; Dudukovic, M. P. Self-Inhibited Rate in Gas-Solid Noncatalytic Reactions. The Shrinking Core Model. Ind. Eng. Chem. Fundam. 1984, 23, 49. Faccio, R. J.; Vidales, A. M.; Zgrablich, G.; Zhdanov, V. P. Connectivity of heterogeneous porous structures and catalytic deactivation. Langmuir 1993, 9, 2499. Giona, M.; Viola, A. The Pressure-Diffusion Equation in the Darcy Regime: Scaling Analysis and Anomalies on Fractals. In Fractals in the Natural and Applied Sciences; Novak, M. M., Ed.; North-Holland: Amsterdam, 1994; pp 165-174. Giona, M.; Adrover, A. Closed-form solution for the reconstruction problem in porous media. AIChE J. 1996, 42, 1407. Giona, M.; Adrover, A.; Giona, A. R. Finite-Difference Methods in the Analysis of Transport Phenomena in Fractal Structures. In Fractals in the Natural and Applied Sciences; Novak, M. M., Ed.; North-Holland: Amsterdam, 1994; pp 153-163. Giona, M.; Adrover, A.; Giona, A. R. Convection-Diffusion Transport in Disordered Structures: Numerical analysis based on the exit-time equation. Chem. Eng. Sci. 1995, 50, 1001. Giona, M.; Schwalm, W. A.; Schwalm, M. K.; Adrover, A. Exact Solution of Linear Transport Equations in Fractal Media I Renormalization Analysis and General Theory. Chem. Eng. Sci. 1996a, 51, 4717. Giona, M.; Schwalm, W. A.; Schwalm, M. K.; Adrover, A. Exact Solution of Linear Transport Equations in Fractal Media II Diffusion and Convection. Chem. Eng. Sci. 1996b, 51, 4731.

Giona, M.; Adrover, A.; Schwalm, W. A.; Schwalm, M. K. Exact Solution of Linear Transport Equations in Fractal Media III Adsorption and Chemical Reaction. Chem. Eng. Sci. 1996c, 51, 5065. Gonzales, A. P.; Pereyra, V. D.; Riccardo, J. L.; Zgrablich, G. The diffusion-controlled annihilation reaction in random adsorptive fields. J. Phys.: Condens. Matter 1994, 6, 1. Gyure, M. F.; Edwards, B. F. Fragmentation of Percolation Clusters at the Percolation Threshold. Phys. Rev. Lett. 1992, 68, 2692. Havlin, S.; Ben-Avraham, D. Diffusion in disordered media. Adv. Phys. 1987, 36, 695. Javier-Cruz, M.; Mayagoitia, V.; Rojas, F. Mechanistic studies of capillary processes in porous media Part II - Construction of porous networks by Monte-Carlo Methods. J. Chem. Soc., Faraday Trans. I 1989, 85, 2079. Joshi, M. Y. A class of stochastic models for porous media, Ph.D. Dissertation, University of Kansas, Lawrence, KS, 1974. Kerstein, A. R.; Bug, A. L. R. Scaling theory of pore growth in a reactive solid. Phys. Rev. B 1986, 34, 1754. Kerstein, A. R.; Edwards, B. F. Percolation model for simulation of char oxidation and fragmentation time-histories. Chem. Eng. Sci. 1987, 42, 1629. Levitz, P. Interfacial geometry of disordered porous media: twopoint correlations and chord length distributions. 4th International Symposium on the Characterization of Porous Solids, Bath, September 15-18, 1996; Book of Abstracts, P4. Mayagoitia, V.; Javier Cruz, M.; Rojas, F. Mechanistic studies of capillary processes in porous media Part I - probabilistic description of porous media. J. Chem. Soc., Faraday Trans. I 1989a, 85, 2071. Mayagoitia, V.; Rojas, F.; Pereyra, V. D.; Zgrablich, G. Mechanistic study of surface processes on adsorbents I. Statistical description of adsorptive surfaces. Surf. Sci. 1989b, 221, 394. Mayagoitia, V.; Rojas, F.; Riccardo, J. L.; Pereyra, V. D.; Zgrablich, G. Dual site-bond description of heterogenoeus surfaces. Phys. Rev. B 1990, 41, 7150. Mohanty, K. K.; Ottino, J. M.; Davis, H. T. Reaction and transport in disordered composite media: Introduction to percolation concepts. Chem. Eng. Sci. 1982, 37, 905. Oran, E. S.; Boris, J. P. Numerical Simulation of Reactive Flows; Elsevier: New York, 1986. Quiblier, J. A. A new three-dimensional modeling technique for studying porous media. J. Colloid Interface Sci. 1984, 98, 84. Ramachadran P. A.; Smith, J. M. Effect of Sintering and Porosity Changes on Rates of Gas-Solid Reactions. Chem. Eng. J. 1977, 14, 137. Reyes, S.; Jensen, K. F. Percolation concepts in the modelling of gas-solid reactions - I. Application to char gasification in the kinetic regime. Chem. Eng. Sci. 1986a, 41, 333. Reyes, S.; Jensen, K. F. Percolation concepts in the modelling of gas-solid reactions - II. Application to char gasification in the diffusion regime. Chem. Eng. Sci, 1986b, 41, 345. Reyes, S.; Jensen, K. F. Percolation concepts in the modelling of gas-solid reactions - III. Application to sulphation of calcinated limestone. Chem. Eng. Sci. 1987, 42, 565. Riccardo, J. L.; Pereira, V.; Rezzano, J. L.; Rodriguez Saa´, D. A.; Zgrablich, G. Effect of adsorptive energy correlation on adsorption isotherms of lattice gases on heterogeneous surfaces. Surf. Sci. 1988, 204, 289. Riccardo, J. L.; Chade, M. A.; Pereyra, V. D.; Zgrablich, G. Adsorption and surface diffusion on generalized heterogeneous surfaces. Langmuir 1992, 8, 1518. Sahimi, M. Applications of Percolation Theory; Taylor & Francis: London, 1994. Sahimi, M.; Tsotsis, T. T. Dynamic Scaling for the Fragmentation of Reactive Porous Media. Phys. Rev. Lett. 1987, 59, 888. Sahimi, M.; Tsotsis, T. T. Statistical modeling of gas-solid reaction with pore volume growth : Kinetic regime. Chem. Eng. Sci. 1988, 43, 113. Salles, J.; Thovert, J. F.; Adler, P. M. Deposition in porous media and clogging. Chem. Eng. Sci. 1993, 48, 2839. Sapag, K.; Bulnes, F.; Riccardo, J. L.; Pereyra, V.; Zgrablich, G. Tracer diffusion on correlated heterogeneous surfaces. Langmuir 1993, 9, 2670. Selim, M. S.; Seagrave, R. C. Solution of Moving-Boundary Transport Problems in Finite Media by Integral Transforms. I Problems with a Plane Moving Boundary. Ind. Eng. Chem. Fundam. 1973a, 12, 1.

Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997 5009 Selim, M. S.; Seagrave, R. C. Solution of Moving-Boundary Transport Problems in Finite Media by Integral Transforms. II Problems with a Cylindrical or Spherical Moving Boundary. Ind. Eng. Chem. Fundam. 1973b, 12, 9. Shah, N.; Ottino, J. M. Effective transport properties of disordered multiphase composites: Application of real-space renormalization group theory. Chem. Eng. Sci. 1986, 41, 283. Sohn, H. Y.; Xia Y.-N. Effects of Nonuniform Distribution of Solid Reactant on Fluid-Solid Reactions. 1. Initially Nonporous Solids. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 386. Sohn, H. Y.; Xia Y.-N. Effects of Nonuniform Distribution of Solid Reactant on Fluid-Solid Reactions. 2. Porous Solids. Ind. Eng. Chem. Res. 1987, 26, 246. Sotirchos, S. V.; Yu, H.-C. Mathematical modelling of gas-solid reactions with solid product. Chem. Eng. Sci. 1985, 40, 2039. Sotirchos, S. V.; Zarkanitis, S. A distributed pore size and length model for porous media reacting with diminishing porosity. Chem. Eng. Sci. 1993, 48, 1993. Taylor, R.; Krishna, R. Multicomponent Mass Transfer; Wiley: New York, 1993. Yortsos, Y. C.; Sharma, M. Application of Percolation Theory to Noncatalytic Gas-Solid Reactions. AIChE J. 1986, 32, 46.

Yu, H.-C.; Sotirchos, S. V. A Generalized Pore Model for Gas-Solid Reactions Exhibiting Pore Closure. AIChE J. 1987, 33, 382. Zgrablich, G. Surface Diffusion of Adsorbates on heterogeneous Substrates. In Equilibria and Dynamics of Gas Adsorption on heterogeneous Solid Surfaces; Rudzinski, W., Steele, W. A., Zgrablich, G., Eds.; Elsevier: Amsterdam, 1997; pp 373-449. Zgrablich, G.; Mayagoitia, V.; Rojas, F.; Bulnes, F.; Gonzales, A. P.; Nazzarro, M.; Pereyra, V.; Ramirez-Pastor, A. J.; Riccardo, J. L.; Sapag, K. Molecular processes on heterogeneous solid surfaces. Langmuir 1996, 12, 129.

Received for review February 24, 1997 Revised manuscript received June 10, 1997 Accepted June 13, 1997X IE970164X

X Abstract published in Advance ACS Abstracts, August 15, 1997.

A Versatile Lattice Simulator for FluidâSolid Noncatalytic Reactions in

Recommend Documents

A Versatile Lattice Simulator for FluidâSolid Noncatalytic Reactions in