Ind. Eng. Chem. Res. 1999, 38, 4985-4992
4985
Shortcut Method for Lumping Diffusion-Reaction Kinetics in Lamellar Systems Alessandra Adrover* and Francesca Giordano Dipartimento di Ingegneria Chimica, Universita´ di Roma “La Sapienza”, via Eudossiana 18, 00184 Roma, Italy
A simple method for evaluating the time evolution of reactant concentrations for diffusionreaction kinetics in lamellar systems is developed. This method provides an efficient lumping strategy for reducing a system of partial differential equations to a low-dimensional system of ordinary differential equations. In the method proposed, the concentration fields are expressed as a linear combination of the solutions calculated in the two limiting cases of infinitely fast and infinitely slow reactions. The solution for the limiting case of infinitely fast reactions is obtained analytically by evaluating the profile of an auxiliary function (corresponding to the difference function in simple bimolecular reactions), accounting for the global stoichiometry, which satisfies a pure diffusion equation. The method is applied to bimolecular reactions and to more complex schemes, such as parallel and competitive-consecutive reactions in lamellar systems, and provides a satisfactory level of agreement with the solution obtained by means of standard numerical methods. 1. Introduction The rate of reactant consumption in homogeneous fluid-phase reactions involving two mutually diffusing fluid species depends on the degree of segregation.1 As a result, the solution of a reaction-diffusion problem in a mixing system requires a model capable of describing mixing and providing some numerical/approximate methods for solving the material balance equations. The special case of the mixing of very viscous fluids is of great interest in many industrial applications, such as polymerization reactions.2 The effects of mixing on the course of such reactions in the presence of diffusion have been analyzed experimentally and numerically for both premixed and unpremixed feeds (see for instance refs 3 and 4 and references therein). Lamellar models, introduced by the Minnesota group,5,6 provide a simplified one-dimensional description of higher-dimensional structures that accounts for the complex striation structure resulting from the stretching and folding dynamics generated by the mixing protocol (for fluids with similar physical properties and no interfacial tension effects). Depending on the characteristic time of reaction, diffusion takes place with the effect of destroying the lamellar structure, as analyzed by Muzzio and Ottino.7 In a diffusion-controlled regime, Sokolov and Blumen8 and Sokolov et al.9 propose a method to solve the material balance equations for bimolecular reactions based on the difference function between reactant concentrations. More recent papers addressing lamellar modeling are due to Clifford and co-workers.10,11 The aim of this work is to propose a simple lumping method for solving diffusion-reaction kinetics in lamellar systems. The concentrations of the reacting species are expressed as a linear combination of the solutions in the two limiting cases of diffusion-controlled and reaction-controlled kinetics. This provides an efficient lump* Corresponding author. Telephone: +39-06-44585892. Fax: +39-06-44585339. E-mail:
[email protected].
ing strategy, as it reduces the numerical solution of a system of partial differential equations (PDEs for short), expressing the material balance for all the reactants to a set of ordinary differential equations (ODEs) for the expansion coefficients. The number of such coefficients is at most equal to the number of reacting species. It is worth noting that the idea of linear interpolation between the two limiting regimes has also been pursued by Baldyga12,13 in the different context of fast homogeneous reactions in tubular reactors in a turbulent regime. This approach is, however, grounded entirely on the introduction of a probability density function of chemically inert species (usually a β function) applied to derive the mathematical expression of the average value of the nonlinear reaction terms in the turbulent mass balances of the reacting species. The article is organized as follows. Section 2 describes the method proposed by considering the simple case of a single-step bimolecular reaction in order to show that a simple macroscopic balance reduces the set of PDEs to a single ODE for a single time-dependent coefficient. Section 3 extends the method to more complex reaction schemes involving more than two reactants, such as parallel and consecutive-competitive reactions, which require only two or three expansion coefficients, respectively. The analysis of complex multilamellar structures is also briefly addressed. The simulation results over a wide range of parameters controlling the reaction conditions reveal a good level of quantitative agreement between the approximated lumping method proposed and the solution of the full set of PDEs. 2. Formulation of the Lumping Procedure In this section, the analysis of diffusion-reaction phenomena in a lamellar system focuses on the simplest case of single-step bimolecular reactions A + bB f products, involving two initially segregated fluids, in order to introduce the basic formulation of the lumping method proposed. Each lamella is initially constituted by two adjacent half striations containing the two
10.1021/ie990353a CCC: $18.00 © 1999 American Chemical Society Published on Web 11/06/1999
4986
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
species A and B. In such a system, diffusion is essentially a one-dimensional problem. The entire lamellar structure may be thought of as the periodic repetition of an elementary lamella of length L. The local material balances for the two reactants in the elementary lamella read as
∂cA ∂2cA ) DA,m 2 - kcAcB ∂t ∂x ∂cB ∂2cB ) DB,m 2 - bkcAcB ∂t ∂x
| |
∂cB ∂x
| |
∂cA ∂x
) x)0
∂cB ∂x
) x)0
x)L
x)L
z(τ,ξ) ) )0
(3)
to be solved with the dimensionless initial and boundary conditions u(0,ξ) ) u0, w(0,ξ) ) w0 and
| |
∂u ∂u ) )0 ∂ξ ξ)0 ∂ξ ξ)1 ∂w ∂w ) )0 ∂ξ ξ)0 ∂ξ ξ)1
w°(τ,ξ) )
b0 2
a0
∞
+
2
ane-λ ∑ n)1
∑
nτ
zne-λ ∑ n)1
2
nτ
cos(λnξ)
b ne
-λ2nτ
(7)
from which it is possible to evaluate the reaction front position δ(τ) at each time (by simply enforcing z(τ,δ(τ)) ) 0) and the reactant concentration profiles
u∞(τ,ξ) ) z+(τ,ξ) )
z + |z| 2
w∞(τ,ξ) ) -z-(τ,ξ) ) -
z - |z| 2
(8)
where z+(τ,ξ) and z-(τ,ξ) are respectively the positive and the negative parts of the z profile. For a generic finite and nonvanishing value of φ2, neither diffusion nor reaction is completely limiting and eq 3 is not susceptible to a closed-form solution. It is in fact possible to introduce a consistent conceptual and practical simplification by approximating the solution of eq 3 as a convex linear combination of the two limiting solutions corresponding to φ2 ) 0 and φ2 f ∞:
(4)
cos(λnξ)
w(τ,ξ) ) β(τ) w°(τ,ξ) + (1 - β(τ))w∞(τ,ξ)
cos(λnξ), λn ) nπ (5)
n)1
where {an, bn; n ) 0, ..., ∞} are constant coefficients directly obtained from the initial conditions.
(9)
As we shall demonstrate, it is only necessary to solve a two-dimensional ODE system in order to obtain R(τ) and β(τ), the other quantities in eq 9 being already available in closed form. Let us refer to the integral quantities U(τ) and W(τ) obtained by integration of the concentration fields over the total length of the lamella. For φ2 ) 0, U° and W° are constant values U°(τ) ) a0/2 and W°(τ) ) b0/2 for all τ g 0, whereas for φ2 f ∞ U∞(τ) and W∞(τ) can be straightforwardly obtained analytically from eq 8. By applying the linear approximation eq 9, it follows that
U(τ) ) R(τ) U°(τ) + (1 - R(τ))U∞(τ) W(τ) ) β(τ) W°(τ) + (1 - β(τ))W∞(τ)
∞
+
2
2
∞
+
u(τ,ξ) ) R(τ) u°(τ,ξ) + (1 - R(τ))u∞(τ,ξ)
Depending on the value of the Thiele modulus (The square Thiele modulus φ2 coincides with the second Damko¨hler number DaII used elsewhere, for example, ref 3.), that is, on the ratio of the characteristic diffusion time tD ) L2/D to the characteristic reaction time tR ) 1/(kcr), different regimes become dominating. For φ2 ) 0, no reaction occurs and eqs 3, corresponding to a pure diffusion problem, can be solved analytically, obtaining the following Fourier-series expansion for the reactant concentrations u°(τ,ξ) and w°(τ,ξ):
u°(τ,ξ) )
z0
(2)
∂u ∂2u ) - φ2uw ∂τ ∂ξ2 ∂w ∂2w ) 2 - bφ2uw ∂τ ∂ξ
(6)
with z(0,ξ) ) u0 - w0/b and ∂z(τ,0)/∂ξ ) ∂z(τ,1)/∂ξ ) 0. Equation 6 can be solved in closed form, thus obtaining the following Fourier-series expansion for the z function:
)0
to be solved with the initial conditions cA(0,x) ) cA0 and cB(0,x) ) cB0. Let us assume equal diffusivities, that is, DA,m ) DB,m ) D. By introducing the dimensionless quantities u ) cA/cr, w ) cB/cr (cr being some reference concentration), ξ ) x/L, and τ ) tD/L2 and the square Thiele modulus φ2 ) kcrL2/D, eqs 1 can be rewritten in dimensionless form
| |
∂z ∂2z ) ∂τ ∂ξ2
(1)
equipped with the boundary conditions for a closed system (due to periodicity)
∂cA ∂x
For φ2 f ∞, that is, for infinitely fast reactions, reaction occurs solely at the interface between the two fluids and the reactants remain completely segregated as the reaction proceeds. In this limiting case, a simple way of evaluating the temporal evolution of the concentration profiles u∞(τ,ξ) and w∞(τ,ξ) and the reaction front position δ(τ) is based on the introduction of the difference function z(τ,ξ) ) u(τ,ξ) - w(τ,ξ)/b, which by definition satisfies the pure diffusion equation
(10)
The expansion coefficients R(τ) and β(τ) entering into eq 9 can be obtained by substituting these expressions in the local material balance equations eq 3 and integrating over the whole system. By substituting eq
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4987
9 into eq 3 and integrating over ξ, we obtain a system of two ODEs for the coefficients R and β:
(U° - U∞)
dU∞ dR ) -(1 - R) dτ dτ
∞
∫01u(τ,ξ) w(τ,ξ) dξ ) -(1 - R) dU dτ
φ2
-
φ2[R(τ) β(τ) H1(τ) + R(τ)(1 - β(τ))H2(τ) + β(τ)(1 - R(τ))H3(τ)] (11) (W° - W∞)
dW∞ dβ ) -(1 - β) dτ dτ
∞
∫01u(τ,ξ) w(τ,ξ) dξ ) -(1 - β) dW dτ
bφ2
-
bφ2[R(τ) β(τ) H1(τ) + R(τ)(1 - β(τ))H2(τ) + β(τ)(1 - R(τ))H3(τ)] with
H1(τ) )
∫0 u°(τ,ξ) w°(τ,ξ) dξ
H2(τ) )
∫01u°(τ,ξ) w∞(τ,ξ) dξ
H3(τ) )
∫
Figure 1. U(τ) versus τ at δ0 ) 0.5 for several values of ζ: (a) ζ ) 0.2; (b) ζ ) 0.5; (c) ζ ) 0.7. Lumped approximation (continuous line) versus numerical solution of eq 3 (]): (A) φ ) 1; (B) φ ) 2; (C) φ ) 10; (D) φ ) 20. The dotted line in part C corresponds to the solution U∞(τ) in the limit for φ f ∞.
1
1 ∞ u (τ,ξ) 0
w°(τ,ξ) dξ
(12)
Special attention should to be paid to the initial condition for eq 11. As a consequence of the hypothesis of initially segregated reactants, the integral ∫10 u(τ,ξ) w(τ,ξ) dξ is vanishing, since A and B do not coexist. Equation 11 at τ ) 0 therefore becomes
1 - R(0) dU∞(0) dR )dτ U° - U∞ dτ
(13)
To avoid the singularity for U° - U∞ at time τ ) 0, it is necessary that R(0) ) 1. Totally analogous considerations apply with regard to the coefficient β(τ), for which we obtain β(0) ) 1. The only set of initial conditions making it possible to avoid an unphysical singularity at τ ) 0 for eq 11 is therefore R(0) ) β(0) ) 1. A further simplification can be obtained by considering the global stoichiometry of the reaction. From a macroscopic material balance it follows that
U(τ) -
W(0) W(τ) W° ) U(0) ) U° ) const b b b
(14)
By substituting eq 10 in this result, it follows that
R(τ) ) β(τ)
(15)
for each τ g 0. A single expansion coefficient, say R(τ), is therefore required for lumping bimolecular reactions. To sum up, the lumping method expressed by eq 9 enables us to determine an approximate solution of generic bimolecular diffusion-reaction kinetics by solving a single ODE to which the original PDE set eq 3 is reduced. The validity of this assumption has been tested by comparison with the solutions of eq 3 obtained by means of a standard numerical method (spectral decomposition and numerical integration of the resulting system of generalized Fourier coefficients). The results are given
in terms of the total quantity U of reactant A versus the dimensionless time τ for several values of φ and for different initial concentrations of reactant B, which we assume to be the limiting species (see Figure 1). The following observations arise from simple physical considerations. First, we expect that the lower the parameter φ, the longer the time required for the limiting reactant to disappear. Let U(0) be fixed. Then, for different initial conditions of the limiting species B given by a parameter ζ
ζ ) (u0δ0 - w0(1 - δ0))/(u0 + w0)
(16)
complete consumption of B corresponds to a different saturation value of reactant A. All the figures show that the expected behavior is accurately reproduced. Quantitative comparison of the solution obtained with the method proposed versus the standard numerical solution indicates a better level of agreement for very low and very high values of φ. This can be easily explained since, in the limit of reaction-controlled or diffusioncontrolled kinetics, the model eq 10 tends to neglect the corresponding complementary terms, namely U∞ and U°, respectively. In an intermediate range of φ, the deviations are more pronounced, as one would expect. In particular, parts A-D of Figure 1 refer to φ ) 1, 2, 10, and 20 and to ζ ) 0.2, 0.5, and 0.7, ζ being the dimensionless difference of total initial quantities for A and B previously defined. To show the difference between U∞(τ) and the solution of the diffusion-reaction scheme for relatively high Thiele moduli, the dotted line in Figure 1C illustrates the corresponding solution in the limit of φ f ∞. A comprehensive quantitative measure of the discrepancy between the approximated lumped solution eq 10 and the solution of the full system of PDE is shown in Figure 2. This figure illustrates the maximum normalized deviation
µmax )
|UL - UN| (1 - ζ)
(17)
(UL ) U(τ) for the lumped approximation, UN ) U(τ)
4988
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
u, w, and y being the dimensionless concentrations for A, B, and C, φ21, φ22 being the Thiele moduli for the two reactions, which we assume to be bimolecular, and the initial conditions being u0, w0, and y0. Zero flux conditions are assumed at the boundaries of the elementary lamella (due to the periodicity at the boundaries of the whole lamellar structure). We further assume that A is initially contained in a single striation and that the ideal mixture B + C is contained in the adjacent one. For the sake of simplicity, equal diffusivities are assumed for A, B, and C in the mixture (reactants + products).
Figure 2. Maximum dimensionless deviation µmax (eq 17) between the lumped approximation and the standard numerical solution versus φ for several values of ζ: (a) ζ ) 0.2; (b) ζ ) 0.5; (c) ζ ) 0.7.
for the standard numerical solution) as a function of the Thiele modulus. Curve a refers to ζ ) 0.2, curve b to ζ ) 0.5, and curve c to ζ ) 0.7. As can be observed, the errors range in all cases from a maximum value for φ = 10, that is, for an intermediate value of the parameter, to a minimum for the extreme values φ ) 1 and φ ) 50. To sum up, the satisfactory agreement between the lumping approximation and the standard numerical solution of eq 3 allows us to conclude that the method proposed constitutes a simple and efficient way to approach a diffusion-reaction problem in lamellar systems undergoing bimolecular reactions. We will show in the next section not only that this method is valid for single reactions but also that satisfactory results can be achieved in the presence of more complex reaction schemes. 3. Extensions to Complex Reaction Schemes and to Multilamellar Systems So far, we have only considered the simplest case of single-step reactions, whereas typical industrial reactions normally undergo extremely complicated mechanisms which can be lumped into the two basic schemes of parallel and consecutive-competitive reactions.1 Schemes such as those presented in this paper are in fact of great importance in the analysis of reactions in the presence of mixing.14 In particular, the conversion of a limiting reactant is often used as a segregation index, that is, as a measure of the degree of mixing achieved by imposing a given velocity field. With this aim, several test reactions, both single and competitive, have been analyzed (see for example refs 15 and 16 and references therein). The approximated lumping method developed in the previous section can be directly extended to such systems with satisfactory results. 3.1. Parallel Reactions. The material balance in dimensionless form for three reactants A, B, and C undergoing a parallel reaction scheme A + bB f P, A + cC f S reads as follows:
∂u ∂2u - φ21uw - φ22uy ) ∂τ ∂ξ2
z(τ,ξ) ) u(τ,ξ) -
w(τ,ξ) y(τ,ξ) b c
∂w ∂ w ) 2 - bφ21uw ∂τ ∂ξ (18)
(19)
and considerations similar to those discussed in the previous section can be enforced. By definition, z(τ,ξ) satisfies the pure diffusion equation eq 6 with the initial condition z(0,ξ) ) u(0,ξ) - w(0,ξ)/b - y(0,ξ)/c and zero flux boundary conditions. Its analytical representation is therefore given by the Fourier-series expansion eq 7. Again, the position of the reaction front δ(τ) can be evaluated as z(τ,δ(τ)) ) 0 . The positive and negative parts of z(τ,ξ) allow us to calculate the concentration u∞(τ,ξ) of reactant A and the total concentration of the mixture B + C:
u∞(τ,ξ) ) z+(τ,ξ) w∞(τ,ξ) y∞(τ,ξ) + ) -z-(τ,ξ) b c
(20)
By simply observing that the ratio w(τ,ξ)/y(τ,ξ) is constant and equal to the initial value w(0,ξ)/y(0,ξ), we are able to distinguish between the two components of the mixture. As regards a generic finite value of φ1 and φ2, the approximated lumping method is readily extended to this reaction scheme by simply adding a new equation, or equivalently a new coefficient γ(τ), to model eq 10 in order to account for the presence of species C. This yields
U(τ) ) R(τ)U° + (1 - R(τ))U∞(τ) W(τ) ) β(τ)W° + (1 - β(τ))W∞(τ) Y(τ) ) γ(τ)Y° + (1 - γ(τ))Y∞(τ)
2
∂y ∂2y ) - cφ22uy ∂τ ∂ξ2
The solutions in the limiting cases of φi ) 0 and φi f ∞, i ) 1 and 2, can be easily computed by means of standard methods. More precisely, in a diffusioncontrolled regime, the auxiliary variable z(τ,ξ) can be introduced
(21)
By integrating eq 18 over the total length of the system and substituting eq 21 in the integral form of eq 18, we obtain a set of ODEs in the coefficients R(τ),
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4989
β(τ), and γ(τ):
dU∞ dR ) -(1 - R) dτ dτ
(U° - U∞) φ21
∫01u(τ,ξ) w(τ,ξ) dξ - φ22∫01u(τ,ξ) y(τ,ξ) dξ ) dU∞ - φ21[R(τ) β(τ) H1(τ) + R(τ) × dτ (1 - β(τ))H2(τ) + β(τ)(1 - R(τ))H3(τ)] -
-(1 - R)
φ22[R(τ) γ(τ) K1(τ) + R(τ)(1 - γ(τ))K2(τ) + γ(τ)(1 - R(τ))K3(τ)] (W° - W∞)
dβ dW∞ ) -(1 - β) dτ dτ
dW∞ bφ 0 u(τ,ξ) w(τ,ξ) dξ ) -(1 - β) dτ 2 bφ1[R(τ) β(τ) H1(τ) + R(τ)(1 - β(τ))H2(τ) + β(τ)(1 - R(τ))H3(τ)]
∫
2
dY∞ dγ ) -(1 - γ) (Y° - Y ) dτ dτ
∞
∫01u(τ,ξ) y(τ,ξ) dξ ) -(1 - γ) dY dτ
-
cφ22[R(τ) γ(τ) K1(τ) + R(τ)(1 - γ(τ))K2(τ) + γ(τ)(1 - R(τ))K3(τ)] (22) with
H1(τ) )
∫01u°(τ,ξ) w°(τ,ξ) dξ, H2(τ) )
H3(τ) )
∫
1 ∞ u (τ,ξ) 0
∫01u°(τ,ξ) w∞(τ,ξ) dξ
w°(τ,ξ) dξ,
∫01u°(τ,ξ) y∞(τ,ξ) dξ, K3(τ) )
∫01u°(τ,ξ) y°(τ,ξ) dξ
∫01u∞(τ,ξ) y°(τ,ξ) dξ
(23)
At time τ ) 0, due to the initial segregation of the reactants, the integrals on the right-hand side of eq 22 vanish, so that, to avoid the singularity at time τ ) 0, the initial condition R(0) ) β(0) ) γ(0) ) 1 must be imposed to solve eq 22. Moreover, it is easy to show that only two coefficients, say R(τ) and β(τ), are needed, since, from the macroscopic material balance it follows that
U(τ) -
ζ)
u0δ0 - (1 - δ0)(w0/b + y0/c) u0 + w0/b + y0/c
(25)
and δ0 are set to 0.2 and 0.5, respectively. In particular, Figure 3 A refers to φ1 ) 20 and φ2 ) 1, 5, 10, and 20, while Figure 3B corresponds to φ1 ) 1 and φ2 ) 1, 5, and 10. As can be observed, the level of agreement between the method proposed and the numerical solution is satisfactory and the maximum normalized deviation is acceptable in both cases. 3.2. Consecutive-Competitive Reactions. Let us now examine the case of three reactants undergoing A + b1/a1B f r1/a1R and R + b2/r2B f S, both bimolecular with constants k1 and k2. We assume that A and B are initially segregated while R is initially (τ ) 0) absent. The local material balance equations are now in dimensionless form
∂u ∂2u ) - φ21uw ∂τ ∂ξ2
K1(τ) ) K2(τ) )
the value 0.5. The other parameters
1
∞
cφ22
Figure 3. Comparison of the approximated lumping solution (solid line) and the numerical solution (]) obtained by the integration of eq 18. (A) φ1 ) 20, δ0 ) 0.5, ζ ) 0.2, (a) φ2 ) 1, (b) φ2 ) 5, (c) φ2 ) 10, (d) φ2 ) 20; (B) φ1 ) 1, δ0 ) 0.5, ζ ) 0.2, (a) φ2 ) 1, (b) φ2 ) 5, (c) φ2 ) 10.
W∞(τ) Y∞(τ) W(τ) Y(τ) ) U∞(τ) ) b c b c W° Y° ) const (24) U° b c
and from the condition w/y ) const, γ(τ) can be expressed as a linear combination of R(τ) and β(τ). Figure 3A,B shows a comparison of the results obtained with the lumping approximation proposed and those calculated by numerically solving the set of material balance equations eq 18. These figures show the time behavior of U(τ) versus τ as the parameters φ1 and φ2 range from low to high values. Reactant A is assumed to be nonlimiting, while w0/y0 is kept fixed at
b2 ∂w ∂2w b1 2 ) 2 - φ1uw - φ22wy ∂τ a1 r2 ∂ξ ∂y ∂2y r1 2 + φ uw - φ22wy ) ∂τ ∂ξ2 a1 1
(26)
with u, w, and y dimensionless concentrations of A, B, and R. The initial and boundary conditions are the same as those in the case analyzed above. The material balance can be easily solved in closed form when no reaction occurs, that is, for φi ) 0; i ) 1 and 2. In the case of infinitely fast reactions (φi f ∞; i ) 1 and 2), from stoichiometry we introduce the z function
z(τ,ξ) ) w(τ,ξ) -
(
)
r1 1 b1 + b2 u(τ,ξ) a1 r2
(27)
Species R is formed at the interface at an infinite rate as a result of the first reaction but immediately consumed by the second one. Its presence cannot therefore be detected and Y∞(τ) ) 0 for all τ g 0. The difference function z(τ,ξ) (eq 27) can be obtained as the solution of eq 6 and the concentrations w∞(τ,ξ) and u∞(τ,ξ) are evaluated from the positive and negative parts of z(τ,ξ), due to the complete segregation deriving from the instantaneous nature of the reactions.
4990
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
Figure 4. Comparison of the approximated lumping solution (solid line) and the numerical solution (]) obtained by the integration of eq 26: (a) U(τ); (b) W(τ); (c) Y(τ); (A) φ1 ) 10, φ2 ) 1, r1 ) 2, r2 ) 1; (B) φ1 ) 5, φ2 ) 1, r1 ) 2, r2 ) 1; (C) φ1 ) 10, φ2 ) 5, r1 ) 2, r2 ) 1; (D) φ1 ) 5, φ2 ) 5, r1 ) 1, r2 ) 2.
For generic fast reactions (φ1 and φ2 being finite and nonvanishing) the approximated lumping method can be applied by simply adding a new equation to eq 10, that is, a new coefficient γ(τ) for R:
Y(τ) ) γ(τ)Y°
(28)
being Y∞(τ) ) 0. By substituting eqs 10 and 28 in the material balance eq 26, previously integrated over the total length of the lamellar system, we obtain a set of three ODEs for the coefficients R(τ), β(τ), and γ(τ), which must solved for the set of initial conditions R(0) ) β(0) ) γ(0) ) 1. Once again this initial condition must be enforced in order to avoid any singularity in the ODE evaluated at τ ) 0, A and B being initially segregated by hypothesis and R initially absent. The analysis is performed by considering φ1, φ2, r1, and r2 as parameters and keeping the initial quantities of A and B fixed. Figure 4A refers to φ1 ) 10 and φ2 ) 1, Figure 4B to φ1 ) 5 and φ2 ) 1, and Figure 4C to φ1 ) 10 and φ2 ) 5. In all cases, r1 ) 2 and r2 ) 1. We observe, as expected, that the higher the ratios φ1/φ2 and r1/r2, the larger the quantity of R formed from the two reactions. The level of agreement with the solution of the full set of material balances eq 26 is acceptable in all cases. As previously observed, maximum deviations correspond to higher values of φ1 - φ2. Figure 4D refers instead to a case in which r1 ) 1 and r2 ) 2, with φ1 ) φ2 ) 5: no significant deviations from the numerical solution of eq 26 can be observed. 3.3. Multilamellar Systems. We shall now extend the analysis and the results of the previous sections to multilamellar structures. Arrays of lamellae are typical features of mixed systems. A large number of striations containing one or more fluid species usually alternate with striations of other species. Experimental and numerical evidence of such structures is widely reported in the literature (see ref 17 for example). The invariant properties characterizing the time evolution of the lamellar structure have been analyzed for two- and three-dimensional flows, such as cavity flows and static mixers. The resulting lamellar model accounts for the
initial disorder (nonuniformity) of the structure by considering striations of random length chosen from a given probability density function (pdf). The scaling properties of lamellar systems have been studied8,9 in the case of bimolecular reactions A + bB f products in a diffusion-controlled regime for strictly periodic and random systems. It is also known7 that the striation thickness distribution collapses onto an invariant curve. The lumping procedure analyzed for periodic lamellar arrays gives satisfactory results also in more complex lamellar systems. To this end, let us consider an array of striations of length chosen at random according to a given pdf. A striation of reactant A is followed by a striation of reactant B, and the two fluid species undergo a single bimolecular reaction A + bB f P occurring at a finite rate, that is, with a finite value of the Thiele modulus φ2. When no reaction occurs (φ ) 0), the lamellar structure is only present initially at time τ ) 0 because of the presence of diffusion. The random nature of the initial structure therefore only affects the initial stage of the process, and the material balance can be solved immediately in closed form. For instantaneous reactions (φ f ∞) in random systems, the time evolution of the total number of lamellae and the displacement of the moving reaction fronts must be considered, since the thinner striations tend to disappear rapidly as a result of the reaction occurring at the interface. The introduction of the z function as z(τ,ξ) ) u(τ,ξ) - w(τ,ξ)/b makes it possible to obtain the position of the interfaces δ(τ) by enforcing the condition z(τ,δ(τ)) ) 0 for each lamella. However, in the application of the lumping approach proposed, the detailed knowledge of the thickness of each lamella is of no use, since the integral quantities U∞(τ) and W∞(τ) appearing in eq 10 can be evaluated from the z function by simply enforcing eq 8. The auxiliary function z still satisfies eq 6. The initial conditions can be easily rewritten by choosing a length for each striation from the given pdf. The solution for z(τ,ξ) is therefore once again a Fourier-series expansion eq 7, and the concentrations of A and B are again given by eq 8. In the general case of fast reactions (φ being finite and nonvanishing), the destruction of the lamellar structure caused by the diffusion is significant, and the effect of a random initial distribution therefore only influences the early stages of the process. Extension of the lumping approximation proposed to this case is thus a straightforward matter, as eq 10 still holds. No further analysis is required, and direct use of the results obtained in section 2 is formally possible. Figure 5 shows the integral quantity U of A versus the dimensionless time τ for some values of φ: namely, φ ) 1, 10, and 30. All the cases refer to a set of 100 initial striations with uniform initial size distribution. Figure 5A refers to stoichiometric initial quantities of A and B, U(0) ) W(0) ) 1, while Figure 5B refers to the case in which B is the limiting reactant: W(0) ) 0.3 and U(0) ) 1. Comparison with the numerical integration of eq 3 reveals that the lumping procedure yields an excellent approximation. Further extensions to multilamellar systems subjected to more complex reactions, such as those discussed previously, are a straightforward consequence of the analysis developed in sections 3.1 and 3.2. 4. Concluding Remarks The analysis developed in this manuscript makes it possible to summarize some characteristic features of
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4991
Figure 5. U(τ) versus τ. Approximated lumping solution (solid line) compared to the numerical solution of eq 3 (]) for a set of 100 striations with initial uniform pdf, for φ ) 1 (curve a), φ ) 10 (curve b), φ ) 30 (curve c): (A) initial stoichiometric mixture U(0) ) W(0) ) 1; (B) nonstoichiometric case with U(0) ) 1 and W(0) ) 0.3.
the approximated lumping method as follows. The method has the advantage of reducing a set of partial differential material balance equations into a system of ordinary differential equations for the coefficients of the linear combination Eq 10. All the terms appearing in this representation, namely U°, W°, U∞, and W∞, are available in closed form from the solution of a pure diffusion equation. Moreover, it has been shown that only one coefficient is needed for single bimolecular reactions. In addition to the conceptual simplicity of the method, simulation results show a satisfactory level of agreement with the solutions of the full set of material balance equations eq 3. The method is fully grounded on simple physical observations and can thus be readily applied to more general and complex reaction schemes, as discussed in section 3. Extension to competitive reactions, such as parallel and consecutive reactions, can be performed with minimum effort by simply adding an additional equation (coefficient) for the third fluid species involved in the reactions. In particular, it has been shown that for such schemes the total number of the resulting ordinary differential equations to be solved is at most equal to the number of reacting species. Finally, the extension to more realistic multilamellar systems has been developed by simply changing the expressions of the initial conditions for each reactant in order to account for the random nature of the structure generated by the mixing protocol. Acknowledgment The authors wish to express their gratitude to M. Giona for invaluable discussions. Notation a, b, c, b1, b2, r1, r2:
stoichiometric coefficients
cA, cB:
reactant concentrations
cA0, cB0:
initial reactant concentrations
cr:
reference concentration
DA,m, DB,m, D:
diffusivities
L:
characteristic length
u, w, y:
dimensionless concentrations
u0, w0, y0:
initial dimensionless reactant concentrations
u°, w°, y°:
dimensionless concentrations for infinitely slow reactions
U°, W°, Y°:
dimensionless amounts for infinitely slow reactions
u∞, w∞, y∞:
dimensionless concentrations for instantaneous reactions
U∞, W∞, Y∞:
dimensionless amounts for instantaneous reactions
z:
concentration difference
z+, z-:
positive and negative parts of the z function
Greek Letters R, β, γ:
coefficients of the linear combination (see eq 9 for example)
δ(τ):
dimensionless position of the reaction interface
δ0:
initial dimensionless position of the reaction interface
ζ:
normalized quantity expressed by eq 16
µmax:
normalized deviation of the lumping solution with respect to the solution of eq 3, as expressed by eq 17
ξ:
dimensionless spatial coordinate
τ:
dimensionless time
φ2, φ21, φ22:
square Thiele moduli
Literature Cited (1) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design; John Wiley & Sons: New York, 1979. (2) Chella, R.; Ottino J. M. Modelling of Rapidly-Mixed FastCrosslinking Exothermic Polymerizations. AIChE J. 1983, 29, 373-382. (3) Chella, R.; Ottino, J. M. Conversion and selectivity modifications due to mixing in unpremixed reactors. Chem. Eng. Sci. 1984, 39, 551-567. (4) Fields, S. D.; Ottino, J. M. Effect of striation thickness distribution on the course of an unpremixed polymerization. Chem. Eng. Sci. 1987, 42, 459-465. (5) Ranz, W. E. Applications of a stretch model to mixing, diffusion and reaction in laminar and turbulent flows. AIChE J. 1979, 25, 41-47. (6) Ottino, J. M.; Ranz, W. E.; Macosko, C. W. A lamellar model for analysis of liquid-liquid mixing. Chem. Eng. Sci. 1979, 34, 877890. (7) Muzzio, F. J.; Ottino, J. M. Evolution of a lamellar system with diffusion and reaction: a scaling approach. Phys. Rev. Lett. 1989, 63, 47-50. (8) Sokolov, I. M.; Blumen, A. Diffusion-controlled reactions in lamellar systems. Phys. Rev. A 1991, 43, 2714-2719. (9) Sokolov, I. M.; Schno¨rer, H.; Blumen, A. Diffusion controlled reaction A + B f 0 in one dimension: the role of particle mobilities and the diffusion-equation approach. Phys. Rev. A 1991, 44, 23882393. (10) Clifford, M. J.; Cox, S. M.; Roberts, E. P. L. Lamellar modeling of reaction, diffusion and mixing in a two-dimensional flow. Chem. Eng. Sci. 1998, 71, 49-56. (11) Clifford, M. J. A Gaussian model for reaction and diffusion in a lamellar structure. Chem. Eng. Sci. 1999, 54, 303-310. (12) Baldyga, J. Turbulent mixer model with application to homogeneous, instantaneous chemical reactions. Chem. Eng. Sci. 1989, 44, 1175-1182. (13) Baldyga, J. A closure model for homogeneous chemical reactions. Chem. Eng. Sci. 1994, 49, 1985-2003.
4992
Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999
(14) Villermaux, J. Mixing effects on complex chemical reactions in a stirred reactor. Rev. Chem. Eng. 1991, 7, 51-108. (15) Fournier, M. C.; Falk, L.; Villermaux, J. A new parallel competing reaction system for assessing micromixing efficiencys experimental approach. Chem. Eng. Sci. 1996, 51, 5053. (16) Mao, K. W.; Toor, H. L. Second order chemical reactions with turbulent mixing. Ind. Eng. Chem. Fundam. 1971, 10, 192197.
(17) Ottino, J. M. The Kinematics of Mixing; Cambridge University Press: New York, 1989.
Received for review May 24, 1999 Revised manuscript received September 13, 1999 Accepted September 16, 1999 IE990353A
