Ind. Eng. Chem. Res. 2007, 46, 3283-3295
3283
Enhancement of Gas Absorption by Sparingly Soluble Fine Particles Reacting Instantaneously with the Dissolved Gas: A Cell Model Vinay A. Juvekar,* Abhijeet A. Joshi, and Rochish Thaokar Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India
In the present work, a cell model is developed that is based on the film theory. In contrast to the existing pseudo-homogeneous models, the present approach retains a discreteness of the particles. The gas-liquid film is divided into cylindrical cells. The particles lie on the axes of the cells, and only those nearest to the interface contribute to the enhancement. The distance of these particles from the interface is exponentially distributed. The boundary element method is used to locate the reaction front. The computed flux is ensembleaveraged over the distribution of the particle locations. The effect of particle diameter and volume fraction of the solid on the enhancement factor is calculated. As the particles become finer, the enhancement factor increases significantly above the estimates of the pseudo-homogeneous models. An approximation to the cell model is also developed and is in good agreement with the rigorous model.
Introduction The problem of gas absorption into slurries that contain fine reacting/catalyst particles is frequently encountered in chemical industry. Examples include absorption of acid gases (e.g., SO2, CO2, H2S, Cl2, HCl, NOx) into a slurry of bases (such as Ca(OH)2, Mg(OH)2, CaCO3), halogenation of latex polymers, catalytic oxidation, hydrogenation, etc. When the particles are finer than the liquid diffusion film, their reactions with the dissolved gas near the gas/liquid interface enhances the rate of diffusion and, consequently, the rate of absorption of the gas. The enhancement is most pronounced when the reaction between the dissolved gas and the solute species is instantaneous. The purpose of the present work is to estimate of the enhancement factor for the instantaneous reaction using an approach that is different from that already existing in the literature. Reaction between a soluble gas-phase species (A) and a liquid-phase species (B) can be represented by the following equation:
A(g) f A(l)
(1a)
A(l) + νB(l) f products
(1b)
The film theory postulates existence of a reaction plane, which divides the liquid diffusion film in two distinct regions, viz., the region 0 < z < λ′ (see Figure 1), where only A exists in the solution, and the region λ′ < z < δ, where only B is present. At the reaction plane, the concentrations of both the species are zero and the diffusion fluxes are in the stoichiometric ratio so that A and B are completely consumed. The concentration profiles of A and B in the film are straight lines. The enhancement factor, E0A, which represents the ratio of actual interfacial flux of A (N*A) relative to the maximum achievable flux in the absence of the reaction, is written as * To whom correspondence should be addressed. Tel.: +91(22) 2576 7236. Fax: +91(22) 2572 6895. E-mail:
[email protected].
Figure 1. Model of Ramachandran and Sharma. Dotted lines represent the reaction front and the concentration profiles in the absence of the particles; solid lines correspond to the presence of the dissolving particles. In the absence of particles, the term [B*] should be replaced by [B0].
E0A )
N* A kL[A*]
(2)
For the instantaneous reaction, the enhancement factor is given by
E0A ) 1 + q
(3)
where q is defined as
q)
DB[B0] νDA[A*]
(4)
Here, DA and DB are, respectively, the diffusivities of A and B in the liquid, [B0] is the concentration of B in the liquid bulk, [A*] is the interfacial concentration of A, and ν is the stoichiometric coefficient of B in the reaction (see eq 1b). For the reaction to be considered to be instantaneous, the following condition should be satisfied:1
10.1021/ie0612212 CCC: $37.00 © 2007 American Chemical Society Published on Web 03/07/2007
3284
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
Ha .1 q
(5)
Ha in the aforementioned equation represents the Hatta number and is defined (for a first-order reaction between A and B) as
Ha )
xDAk2[B0]
(6)
kL
where k2 is the rate constant of the reaction, and kL is the mass transfer coefficient of A across the gas/liquid interface. According to the film theory kL is given by the ratio DA/δ. We can rewrite the ratio on the left-hand side of eq 5 as
Ha ) q
δ (DAq/xDAk2[B0])
)
δ lr
(7)
where lr can be interpreted as the length of the reaction zone to which the reaction between A and B is confined. Therefore, eq 5 implies that the reaction can be considered instantaneous only when the reaction zone is much thinner than the film. Based on the analysis presented by van Krevelen and Hoftijzer2 in the form of plots of the enhancement factor E0A versus the Hatta number Ha, with q as the parameter, the instantaneous reaction limit (given by eq 5) is determined to be almost attained when the ratio Ha/q is ∼30, which seems to be a reasonable number, given the interpretation presented in eq 7. When solute B is supplied by the dissolution of sparingly soluble fine particles present in the liquid, the reaction is represented by the following equations:
B(s) f B(l)
(8a)
A(g) f A(l)
(8b)
A(l) + νB(l) f products
(8c)
This problem was analyzed for the first time by Ramachandran and Sharma,3 who used the film theory. The dissolution of the particles is assumed to be sufficiently fast to maintain the bulk solution at saturation, with respect to B. If the particles are finer than the thickness of the film, their dissolution in the film modifies the concentration profiles of both A and B (see Figure 1), according to the following equations:
DA
d2[A] dz
2
)
ksap[B*] ν
(for 0 < z < λ)
(9)
and
DB
d2[B] dz2
) ksap([B*] - [B])
(for λ < z < δ) (10)
terms ks and ap represent the solid-liquid mass-transfer coefficient and the specific solid-liquid interfacial area, respectively. The terms on the right-hand side (RHS) of eqs 9 and 10 represent the local rate of dissolution of the particles, where the driving force for dissolution is the difference between the saturation concentration and the local concentration. The aforementioned equations yield the following expressions for the interfacial flux of A:
N*A ) -DA
DA[A*] ksap[B*] λ d[A] |z)0 ) + dz λ ν 2
()
(
) ()
δ-λ 1 1 ) [B*](xDBksap) coth D ka + k a [B*]λ ν DB x B s p ν s p (13) The aforementioned set of equations is first solved for the unknown λ, and N*A can then be obtained from one of the equations, after substituting the value of λ. The enhancement in absorption due to the presence of the particles is given by
EA )
N*A
( )
We can define the relative enhancement factor as
ηA )
[A] ) [B] ) 0 [B] ) [B*]
(at z ) 0) (at z ) λ) (at z ) δ)
EA
(15)
E0A
where E0A is as defined by eq 3. The relative enhancement factor is a measure of the effectiveness of the particles in enhancing the rate of gas absorption. It is always >1, and the larger its magnitude is, the more effective the particles. A somewhat ambiguous quantity in the aforementioned model is ks, which represents the solid-liquid mass-transfer coefficient defined inside a gas-liquid film. This, in the parlance of the film theory, implies diffusion film within another diffusion film! Pal et al.4 recommended the value of ks, which yields a Sherwood number of unity (Sh ) 1), based on the particle radius Rp (ShRp ) ksRp/DB ) 1). This is the Sherwood number that is attained if a spherical particle is placed in an infinite stagnant liquid. This assumption is strictly valid in the limit when the ratio of film thickness to particle diameter in infinite. It is also requires that the interparticle spacing be much larger than the particle diameter, so that neighboring particles do not interact. These conditions are satisfied reasonably well only by very fine particles at very low particle loadings. Using the relation ks ) DB/Rp and writing ap in terms of the volume fraction of the solid in the slurry (φp),
The following boundary conditions prevail:
[A] ) [A*]
(14)
DA [A*] δ
ap )
(11a)
3φp Rp
(16)
(11b) (11c)
we can recast eq 13 in the following dimensionless form:
(1 qδ+ˆ q)(λˆ1q + 3λ2ˆ) qδˆ )( { 3φ coth[(δˆ - λˆ )x3φ ] + 3λˆ φ } 1 + q) x
ηA )
and
DA
DB d[B] d[A] |z)λ ) | dz ν dz z)λ
(12)
Here, [B*] represents the saturation concentration of B. The
p
where
p
p
(17)
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3285
δ)
δ Rp
(18a)
λˆ )
λ Rp
(18b)
and
q)
DB[B*]
(18c)
νDA[A*]
Note that the definition of q given above is obtained from eq 4 by replacing [B0] by [B*], because the bulk solution is saturated, with respect to B. Uchida et al.5 modified Ramachandran and Sharma’s model by accounting for the enhancement in the rate of dissolution of the particles in the region of 0 < z < λ, because of the presence of reactant A. Under this condition, all the equations, except eq 9, remain the same as those given by Ramachandran and Sharma. Equation 9 is modified to
DA
d2[A] dz
2
)
(
)
ksap[B*] ν[A] DA 1+ ν [B*] DB
(for 0 < z < λ) (19)
The term in the parentheses represents the enhancement caused by the instantaneous reaction between A and B near the particle/ liquid interface. The expression for the relative enhancement factor by the model of Uchida et al. is obtained from the following equations:
ηA )
)
( (
){( ){
qδˆ x3φp 1+q
1+
}
1 1 coth(λˆ x3φp) q sinh(λˆ x3φp)
)
qδˆ x3φp coth[(δˆ - λˆ )x3φp] + 1+q
(2 + q1)[coth(λˆ x3φ ) - sinh(λˆ1x3φ )]} (20) p
p
The model of Uchida et al. predicts a higher value of the relative enhancement factor than that based on the Ramachandran and Sharma model. Scala6 incorporated the effect of particle-particle interactions in the film by placing the particle in spherical cells of the fluid. The size of each cell is computed by assigning an equal volume of fluid per particle. Outside the cell, the concentration of A (or B) is uniform and corresponds to the local concentration. The diffusion of A and B inside the cell is considered for determining the rate of solid dissolution. This rate is used in the model of Uchida et al.5 to compute the rate of gas absorption. In effect, this model uses a solid-liquid diffusion film thickness that is equal to one-half of the mean interparticle spacing. The relative enhancement factor by the model can be obtained from the model of Uchida et al. (see eq 19) by replacing φp with φp/(φp-1/3 - 1). The Scala model predicts a lower enhancement factor value than Uchida’s model for φp < 0.125. The reverse is true for φp > 0.125. Uchida et al.,7 Sada et al.,8 and Scala6 extended the theory to the case where the bulk is not saturated, with respect to B. Uchida et al.9 used penetration theory to analyze this problem and obtained an analytical solution for equal diffusivities of the species. Mehra et al.10 analyzed the problem of absorption in a quiescent slurry, where the reaction plane moves with time.
Scala6 also accounted for the presence of gas-film resistance in computing the enhancement factor. All the models previously described can be classified as pseudo-homogeneous, i.e., they assume the solid phase to be microhomogeneously mixed with the liquid in the film, so that the discrete nature of the particles is completely lost. Although the assumption of homogeneity is reasonable in the region a few particle diameters away from the gas/liquid interface, very near the interface, the geometry of the particle, relative to that of the interface, is important for determining the enhancement. For example, consider a spherical particle in the vicinity of a planar interface. Different parts on the surface of the sphere will be located at different distances from the interface. Component A, which is diffusing from the gas/liquid interface into the liquid, will meet component B, which is diffusing outward from the particle surface, at different locations. Thus, we have a curved reaction front in the vicinity of the particle. Because of the fact that particles that are very close to the interface are the ones that distort the reaction front the most, they contribute most to the enhancement. Therefore, the details (such as shape of the particle, its orientation, and the distance of the particle from the interface) are important, in addition to particle size and particle loading. The inherent limitation of the pseudo-homogeneous models is that they cannot account for these details. Moreover, the pseudo-homogeneous models indirectly imply that particles cannot approach the interface closer than the thickness of the solid-liquid film. Thus, in the model that uses ShRp ) 1, the distance of the closest approach equals Rp. According to Scala’s model, the distance of closet approach decreases as the particle loading increases. No such limitation exists in reality and particles can approach arbitrarily close to the gas/liquid interface. These limitations of the existing models can be overcome using an approach that treat particles as distinct entities placed in the gas-liquid diffusion film and acting as sources of B. In the present work, we have developed a cell model to achieve this goal. The film is divided into cylindrical cells of equal volume. We consider spherical particles that are placed on the axes of these cells. The statistical distribution of the distances of the particles from the interface has been postulated. The reaction front is a now a curved surface, rather than a plane. We use the boundary element method to locate this front. The interfacial flux of A and the enhancement factor is then calculated as a function of the distance of the particle from the interface. Based on the distribution of the distances, the statistical average value of the enhancement factor is computed. The effects of the particle diameter, particle loading, ratio of the interfacial concentration of the reactants etc, have been evaluated. Finally, the enhancement factor is compared with the Uchida et al.5 model. The present model can be adapted to penetration theory. Also, a variety of factorssthe gas-film resistance, shape of the particle, change in the particle size due to dissolution, incomplete saturation of the bulk with respect to B, etc.scan also be incorporated into the model. This model would serve as a benchmark to test the accuracy of pseudo-homogeneous models. The Cell Model In the cell model, cells are used to account for the interactions among the particles. When two particles are near each other, each particle prevents spreading of the solute dissolving from the other particle. If all particles are of equal size, and if the interactions among them are symmetric in the mean sense, it is possible to view the particles as being enclosed in equal-volume
3286
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
the interval 0 and Zp does not, in any way, influence its chance of being found in some interval ∆Zp beyond Zp, and (ii) the probability of finding a particle in a short interval ∆Zp varies linearly with ∆Zp. Property (i) is obeyed by the Poisson processes, where the number of events is dependent on the size of the time interval in which they occur, rather than on the location of the time interval on the time axis. In the present problem, the time interval translates to the interval of distance on the cell axis and the number of events corresponds to the number of particles encountered in that interval. Because particles have no preference for a particular location of the interval, property (i) is obeyed in the present case. Property (ii) follows from the general principle that any functional dependence between the probability and ∆Zp can be linearized and this linear approximation becomes increasingly accurate as ∆Zp becomes smaller. It immediately follows from eq 22 that Figure 2. Layout of cells and particles: (a) particles lying on the common axis of cells and (b) two extreme locations of a particle in a cell.
(and equal-shape) cells whose boundaries are reflecting. The total volume of the cells equals that of the liquid. This considerably simplifies the analysis, because, now, only one cell must be examined. The present cell model is similar to that presented by Karve and Juvekar11 for estimating the enhancement by fine catalytic particles at the gas/liquid interface. As the first step, volume of the entire liquid is compartmentalized into equally sized cylindrical cells that have a height-to-diameter ratio of unity. The radius of the cell is denoted by Rc. On average, there is one particle per cell. If the radius of the particle is Rp, and the volume fraction of particles in the slurry is φp, then the radius of the cell can be calculated as
( )
Rc 2 ) Rp 3φp
1/3
P{particle not found within the interval Zp and Zp + ∆Zp|particle not found within the interval (0,Zp)} ) 1 - λ∆Zp (23) Using the relationship between the conditional probability and the joint probability, we can recast eq 23 as
P{particle not found within the interval (0,Zp + ∆Zp)}
) P{particle not found within the interval (0,Zp)} 1 - λ∆Zp (24)
Denoting the probability P{particle not found within the interval (0,Zp)} with the symbol G(Zp), we can rewrite eq 24 as
G(Zp + ∆Zp) G(Zp)
(21)
We assume that the particles lie on the axes of the cells. However, the location of particles on the axis is random. Some cells may be empty, and some may contain more than one particle. Consider now a stack of cells that have a common axis that is normal to the interface. The interparticle spacing is denoted by d, as shown in Figure 2. The particles are randomly distributed on this axis, so that, on average, the interparticle spacing is 2(Rc - Rp). The particle that is nearest to the interface will contribute the most to the enhancement of the absorption of A. Because this particle is encountered first, as we move from the interface into the liquid, we call it the first particle. We are interested in finding the probability distribution of the distance Zp of the leading point of the first particle from the interface. This is done as follows. Consider an observer that starts from the interface and moves along the common axis of the stack of cells, in search of the first particle. If he does not encounter the particle up to a distance Zp from the interface, then a probability P that he would encounter it within a short interval ∆Zp is assumed to be equal to λ∆Zp, where λ is a positive constant that is independent of Zp. Mathematically, we can write
P{particle found within the interval Zp and Zp + ∆Zp|particle not found within the interval (0,Zp)} ) λ∆Zp (22) There are two implications of the aforementioned equation, viz., (i) the fact that the particle has not been encountered within
) 1 - λ∆Zp
(25)
Rearranging the aforementioned equation and applying a limit of ∆Zp f 0, we obtain the following differential equation for G(Zp):
dG(Zp) ) -λG(Zp) dZp
(26)
We solve the equation using the boundary condition G(0) ) 1. This boundary condition implies that the particle is not adsorbed at the gas/liquid interface and, hence, its probability of being located exactly at Zp ) 0 is zero. The solution of eq 26 is
G(Zp) ) exp(-λZp)
(27)
The probability density of location of the first particle, p(Zp), can be obtained from G(Zp) as
p(Zp) ) -
dG(Zp) ) λ exp(- λZp) dZp
(28)
The ensemble-averaged distance of the first particle can now be obtained as
Z hp )
∫0∞ Zpλ exp(-λZp) dZp ) λ1
(29)
It is reasonable to assume that the average distance of the first particle from the interface is the same as the average interparticle distance, i.e., 2(Rc - Rp). Equation 29 then yields λ ) [2(Rc Rp)]-1 and, hence,
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3287
p(Zp) )
[
Zp 1 exp 2(Rc - Rp) 2(Rc - Rp)
]
(30)
We assume that only the first particle contributes to the enhancement of the absorption of A. The particle that has the greatest contribution to the enhancement factor, after the first particle, is the second particle. We expect the contribution due to the second particle to be small, for two reasons. First, it lies beyond the distance Zp + 2Rp from the interface. Hence, it cannot approach the interface as closely as the first particle does. Second, the flux of A, approaching the second particle, is screened by the first particle. We have shown later, through simulations, that even if we neglect the screening effect and consider only the distance effect, the contribution of the second particle is small. We divide a diffusion film that has a thickness δ into equally sized cylindrical cells. Each cylinder has a radius Rc (obtained from eq 21) and a height δ (the film thickness). The particles are placed on the axes of the cells. Depending on the height δ, each newly defined cell may contain one or more of the previously defined cells (of height 2Rc). However, this information is unimportant to us, because we are only interested in the position of the first particle in the new (composite) cell, which is still given by eq 30. If the interfacial flux of A in a cell in which the leading point of the particle is located at a distance Zp from the interface, denoted as N*A(Zp), then the representative flux of A is the ensemble-averaged flux over the distribution of Zp, i.e.,
N*A )
∫0
δ-2Rp
N*A(Zp)p(Zp) dZp + ∞ ∫δ-2R
p
DA[A*](1 + q) p(Zp) dZp (31) δ
The first integral corresponds to the situations when the particle is present in the cell, and the second integral corresponds to the absence of the particle in the cell (where the flux is contributed solely by diffusion from the bulk). The correct upper limit of the first integral (and the lower limit of the second) should be Zp ) δ (the particle is completely out of the film). However, we have used the limit corresponding to a situation when the lowermost point on the particle touches the boundary of the film (see Figure 2b). Thus, we have ignored the contribution to enhancement from the particles that are partly in the film and partly in the bulk. This is done to simplify the computational algorithm. Because these particles are far away from the interface, their contribution to the total enhancement is small, and, hence, we expect the error introduced by this simplification to be small. Figure 3 shows the schematic diagram of a cell. The axis of the cell forms the z-axis. The top plane surface of the cylinder, which is also the gas/liquid interface, corresponds to z ) 0. The bottom plane (z ) δ) corresponds to the lower boundary of the diffusion film. The radial coordinate spans from r ) 0 (axis) to r ) Rc (the curved boundary of the cell). The azimuthal angle spans from φ ) 0° to φ ) 2π. The center of the particle lies on the axis of the cell at z ) Zp + Rp, so that its apex D lies at z ) Zp. The projection of the reaction front is shown by the curve AB. The front itself is the curved surface obtained by revolving AB about the cell axis. The front divides the cell in two regions. In region 1, which is enclosed by the boundary CABIC, only the reactant A is present, whereas in region 2, which is enclosed by ADEFGHBA, only B is present. On the front AB, both A and B react instantaneously, so that their concentrations fall to zero.
Figure 3. Details of the cell and the boundary conditions. Cross-hatched region (region 1) contains only component A, and the remainder (region 2) contains only component B.
The equation governing the diffusion of A in region 1 is the Laplace equation,
∇2[A] ) 0
(32)
where
∇2 ≡
1 ∂ ∂ ∂ r + 2 r ∂r ∂r ∂z
( )
Note the absence of an azimuthal coordinate φ, which is due to the axial symmetry. The boundary conditions are
on the boundary CI: on the boundary AB:
[A] ) [A*] [A] ) 0
on the boundaries CA and IB:
(33a) (33b)
∂[A] )0 ∂r (33c)
The Laplace equation also governs the diffusion of B in region 2:
∇2[B] ) 0
(34)
The following boundary conditions must be imposed:
on the boundaries GH and DEF: on the boundary AB:
[B] ) [B*]
[B] ) 0
on the boundaries AD, FG, and HB:
(35a) (35b)
∂[B] )0 ∂r
(35c)
The boundary condition on GH is based on the assumption that the bulk (z g δ) is saturated, with respect to solute B. In the boundary conditions given by eqs 33c and 35c, we assume that the fluxes of A and B across the curved boundary of the cell are zero. This is the result of the basic assumption of the cell model that the curved boundary of the cell, which expresses the interaction among the particles, is a reflecting boundary. This assumption tacitly implies that the interaction among the neighboring particles is symmetric. In the strict sense, symmetry requires that all the neighboring first particles are located the same distance from the interface as that of the particle under consideration. This is not true, because the distances of the neighboring first particles from the interface
3288
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
are randomly distributed. However, we can justify this assumption in the ensemble-averaged sense. Because there is no spatial variation of particle concentration along the interface, on average, all the cells are equivalent. Hence, no cell can preferentially leak, either species A or B, into the neighboring cell through any part of the curved boundary. As a result, although the actual value of N*A(Zp) would fluctuate from cell to cell, about the value computed by the present procedure, we expect these fluctuations to cancel during the ensemble averaging by eq 31 to yield negligible residual deviation. The aforementioned set of equations can be simultaneously solved to obtain the concentration profiles of A and B in regions 1 and 2, respectively, provided that the location of the reaction front AB is known a priori. However, because of the fact that we do not have this knowledge, we must track the reaction front. The information to do so is provided by the flux balance. On the boundary AB, the fluxes of A and B should be opposite in direction and in stoichiometric proportion, i.e.,
DA
DB ∂[B] ∂[A] | ) | ∂nj1 AB ν ∂nj2 AB
(36)
N*A(Zp) ) -
2
πRc
∫0R
c
∂[A] | (2πr) dr ∂z z)0
N*A )
1 2(Rc - Rp)
∫0δ-2R
p
N*A(Zp)e-Zp/[2(Rc-Rp)] dZp +
DA[A*] (1 + q)e-(δ-2Rp)/[2(Rc-Rp)] (38) δ
δN* A(Zp)
A ˆ )1
DA[A*](1 + q)
The ensemble average relative enhancement factor ηA (defined as ηA ) δN*A/{DA[A*](1 + q)}) can be computed using the following equation, which is derived from eq 38:
ηA )
1 2(Rc - Rp)
∫0δ-2R
p
ηA(Zp)e-Zp/[2(Rc-Rp)] dZp + e-(δ-2Rp)/[2(Rc-Rp)] (40)
Method of Solution Equations 32 and 34 and the corresponding boundary conditions are first converted to dimensionless form, using the following transformations:
[B] [B*]
(41b)
rˆ )
r Rp
(41c)
zˆ )
z Rp
(41d)
Rc )
Rc Rp
(41e)
δˆ )
δ Rp
(41f)
Zˆ p )
Zp Rp
(41g)
(42)
at z ) 0, 0 < rˆ < Rˆ c
(43)
A ˆ )0
at zˆ ) zˆλ (rˆ), 0 e rˆ e Rˆ c
(44)
∂A ˆ )0 ∂rˆ
at rˆ ) 0, 0 < zˆ < zˆλ(0)
(45)
at rˆ ) Rˆ c, 0 < zˆ < zˆλ (Rˆ c)
(46)
∂A ˆ )0 ∂rˆ
In the aforementioned equations, zˆλ(rˆ) is the dimensionless z-coordinate of the reaction front with zˆλ(0) and zˆλ(Rˆ c) being the corresponding values on the cell axis (rˆ ) 0) and the curved boundary of the cell (rˆ ) Rˆ c). In region 2:
ˆ ∂B ˆ 1 ∂ ∂B rˆ + 2)0 rˆ ∂rˆ ∂rˆ ∂zˆ
( )
B ˆ )1 B ˆ )1
(39)
B ˆ)
( )
The relative local enhancement factor in the region of 0 e Zp e δ - 2Rp is obtained from the following expression:
ηA(Zp) )
(41a)
ˆ ∂A ˆ 1 ∂ ∂A rˆ + 2)0 rˆ ∂rˆ ∂rˆ ∂zˆ
(37)
The aforementioned flux is a function of Zp, the particle location. This flux is ensemble-averaged, using eq 31, which can be rewritten using eq 30 as
[A] [A*]
The dimensionless forms of the equations and the boundary conditions are as follows: In region 1:
nj1 and nj2 are the unit normals to AB; nj1 points into region 2 from region 1 and nj2 points into region 1 from region 2. Note that nj2 ) -nj1. The area-averaged flux of A crossing the interface at z ) 0 is
DA
A ˆ )
∂B ˆ )0 ∂rˆ
(47)
at zˆ ) δˆ , 0 < rˆ < Rˆ c
(48)
at rˆ ) rˆp (zˆ), Zˆ p e zˆ e Zˆ p + 2
(49)
at rˆ ) 0, zˆλ(0) < zˆ < Zˆ p, and Zˆ p + 2 < zˆ < δˆ (50) ∂B ˆ )0 ∂rˆ
at rˆ ) Rˆ c, zˆλ(Rˆ c) < zˆ < δˆ
(51)
In the aforementioned equations, rˆp(zˆ) is the dimensionless rˆcoordinate of a point on the surface of the particle, written as a function of the zˆ-coordinate. Equation 36 for the reaction front is transformed to
∂B ˆ| ∂A ˆ| | ) -q | ∂nˆ 2|zˆλ ∂nˆ 2|zˆλ
(52)
The procedure to locate the front is iterative. An initial location is guessed. Equations 42 and 47 are then solved and the
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3289
mismatch between the two sides of eq 52 is computed at various points on the front. Depending on the magnitude of the mismatch at a point, the point is appropriately moved along the direction of the normal. After the current location of the front is determined, the procedure is repeated until the mismatches are reduced to within the acceptable level. We have used the boundary element method (BEM) to solve these equations. The Laplacian nature of the equations makes the use of BEM particularly convenient. The method was selected because of the now well-appreciated advantages, such as reduction of the computation domain and more-accurate solutions over other discretization methods. We first convert the equations to the boundary integral equation (BIE) form.12 The BIEs of the axisymmetric Laplace equation for regions 1 and 2 are
A ˆ (p) + 2π
∫Γ
1
K1(p, Q)A ˆ (Q)rQ dΓQ )
K1(p,Q) )
[
n1r 1 π K m, 2 2 2 2π C 2π rQ
( ) 2(z - z ) F π π n E(m, ) n E(m, ) (60) 2 D 2] 2π r D p
1z
Q
where
D ) [(rp2 - rQ2) + (zp - zQ)2]
(61a)
F ) rp2 - rQ2 + (zp - zQ)2
(61b)
and E(m,π/2) is the elliptical integral of the second type:
( π2) ) ∫
π/2
E m,
0
∫
∫Γ
2
K1(p, Q)B ˆ (Q)rQ dΓQ ) ∂B ˆ (Q) 2π Γ K2(p, Q) r dΓ (54) 2 ∂n Q Q
∫
Γ1 and Γ2 are the boundaries of regions 1 and 2, respectively. The index “p” corresponds to a point, in the domain or on the boundary, at which the value of A ˆ or B ˆ is desired and Q is a variable point that is moved along the boundary. The term rQ represents the r-coordinate of point Q. The integrals in these equations are the line integrals (boundary integrals), to be performed along the boundaryΓ1. K1(p,Q) and K2(p,Q) are the kernels of the respective boundary integrals. K2(p,Q) is defined as
K2(p,Q) )
1 π K m, 2 2 2π C
( )
(55)
where C is defined as
C ) x(rp + rQ)2 + (zp - zQ)2
(56)
ri and zi are the r- and z-coordinates of point i (where i ) p or Q) and K(m,π/2) is the elliptic integral of the first type, which is defined as
( π2) ) ∫
K m,
dR
π/2
0
x1 - m2 sin2 R
m)
4π2C
K1|rp)0 ) K2|rp)0 )
∂K2(p,Q) ∂n1
n1r )
∂r ∂n1
(63a)
n1z )
∂z ∂n1
(63b)
rQn1r - (zp - zQ)nz 2π[rQ2 + (zp - zQ)2]3/2 2
2π[rQ
-1 + (zp - zQ)2]1/2
(64) (65)
The numerical procedure that is used to solve the boundary integrals in eqs 53 and 54 involves discretization of the boundary into several connected elements. We use two nodded, equalarc-length boundary elements. The points, which are internal to the element, are expressed in terms of the end nodes, using suitable interpolation polynomials. For straight boundaries, one of the coordinates remains constant and no interpolation is required. For the curved boundaries (AB and DEF), we use cubic spline interpolation.13 The concentrations A ˆ and B ˆ on the elements are linearly interpolated. Thus, the dimensionless concentration of A at a point in an element, the nodes of which are j and j + 1, is given by
A ˆ (r,z) ) Rrˆ + βzˆ
(57)
(66)
where
Ajzˆj+1 - Aj+1zˆj rˆjzˆj+1 - rˆj+1zˆj
(67a)
- Ajrˆj+1 + Aj+1rˆj rˆjzˆj+1 - rˆj+1zˆj
(67b)
R) (58)
K1(p,Q) is the normal derivative of K2(p,Q) on the boundary,
K1(p,Q) )
(62)
and the expressions for the kernels given by eq 55 and 59 are valid for the points that do not fall on the axis of the cell. For the points on the axis (i.e., on the line r ) 0), the following expressions must be used for the kernels:
with m given by the following expression:
2xrprQ
x1 - m2 sin2 R dR
The components n1r and n1z of the normal nj1 are given by
∂A ˆ (Q) 2π Γ K2(p, Q) r dΓ (53) 1 ∂n Q Q B ˆ (p) + 2π
Q
1r
2
(59)
where n1 is the outward normal to region 1. We can explicitly obtain K1(p,Q) as
and
β)
A similar equation can be written for B ˆ. The integrals in eqs 53 and 54 are evaluated using a 16-point Gaussian quadrature formula. The integrals show a logarithmic singularity14 when the point p approaches Q (i.e., m f 1), because, under these conditions, K(m,π/2) f
3290
Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007
Figure 4. Effect of the cell radius on the profile of the reaction front: (a) Rˆ c ) 2, (b) Rˆ c ) 4, and (c) Rˆ c ) 8. In all panels, q ) 1 , δˆ ) 5, and Zˆ p ) 0.774.
ln(4/x1-m2). The singular integrals are evaluated by subtracting the singular part and evaluating the nonsingular integrals by ordinary Gaussian quadrature. The singular part is then analytically integrated, assuming straight elements. The numbers of elements used were as follows: 40 on arc AB, 40 on arc DEF, 40 each on the straight boundaries GH and CI, 10 on boundary BI, and 20 on the boundary HB. The convergence criterion used was that the highest magnitude of the difference in the fluxes of A and B on the boundary AB should be
