Ind. Eng. Chem. Res. 1995,34, 794-805
794
Reaction Fronts in a Porous Medium. Approximation Techniques versus Numerical Solution Fernando Escobedo and Hendrik J. ViUoen* Department of Chemical Engineering, University of Nebraska-Lincoln, Lincoln, Nebraska 68588-0126
The flame sheet approximation (FS) and a novel polynomial approximation technique (PA) are compared in terms of their capability to describe reaction fronts of highly exothermic reactions in a porous medium. A one-phase model and a two-phase model of a system with adiabatic walls and a radiant output (to approximate the case of a porous radiant burner) are included in the analysis. By matching the reaction zone solution found by either the FS or PA method with the solutions of the nonreacting zones, the temperature, conversion, and position of the reaction zone were determined. Numerical solutions for catalytic and noncatalytic oxidation reactions were used to compare the predictions of both approaches. It was found that although both techniques yielded good approximations to the solutions, the PA technique proved t o be more accurate, producing results with 3.5% of the numerical results. Both methods can find useful application in the analysis of this class of problems. 1. Introduction
Exothermic reactions are commonly carried out in catalytic reactors in many practical applications. In some cases, (i.e., oxidation of methane and carbon monoxide) the reaction is carried out in the ignited state. The combustion of a gas fuel inside a porous media, as found in radiant porous burners (PRB), is one of the most promising novel combustion concepts (Viskanta, 1991). These burners stabilize a premixed flame inside a noncatalytic porous medium. The heat of combustion released in the gas phase heats the porous matrix by convection which in turn heats the load by radiation, increasing the overall heat transfer efficiency (Viskanta, 1991). In all these cases the high activation energy of the reaction limits the reaction zone to a rather thin section of the bed. By matching the (approximate) solution inside the reaction zone with the solutions outside the reaction zone, analytical expressions may be obtained for the “flame” position, temperature, and conversion in terms of other parameters of the system. The reaction zone has been successfully described by the “flame sheet” approximation in various applications. Gatica et al. (1987) solved numerically the problem of traveling waves in nonadiabatic packed bed reactors. On the basis of the flame sheet approximation, they also provided expressions for the velocity of the traveling waves. This method has been successfully applied to other combustion problems (cf. Matkowsky and Sivashinsky (1979)and Norbury and Stuart (1989)). In all these applications, the system always consisted of an infinitely long reactor. A finite axial geometry has only recently been analyzed (Escobedo and Viljoen, 1993). Various authors have modeled PRB’s with a different degree of sophistication (Takeno et al.,1981; Yoshizawa et al., 1988; Singh et al., 1991; Sathe et al., 1990; etc.). A numerical approach has been always used to solve the model. In all the cases the position of the flame was assumed to be known a priori and the burning rate (fuel flux) was computed. An analytical treatment of the most rigorous expressions for the radiation heat transfer and for the mass
* Author to whom correspondence should
be addressed.
0888-588519512634-0794$09.00~0
and energy balances is still missing, and researchers relied on numerical methods to analyze the problem. The use of simplified models leads to more tractable problems which still capture the essential features of the physical phenomena. It also simplifies optimization, parameter sensitivity, or parameter identification problems. The simplifications which have been used include the assumption of thermal equilibrium between phases, constant solid phase temperature (Takeno et al., 19811, the use of the delta function to describe the reaction rate (Singh et al., 19911, approximation of the radiation heat transfer by the two-flux model (Singh et al., 19911, or the diffusion model (Lawson and Norbury, 1985). The analysis which was presented by the authors in a previous paper (Escobedo and Viljoen, 1993) is extended. A matched asymptotics approach is also followed to solve a one-phase (pseudohomogeneous)and a two-phase model. An alternative method to describe the reaction zone is also proposed. Polynomial approximations are used to represent the temperature profiles in that zone. This approach serves as a transition between the analytical attempts to solve the models approximately and purely numerical techniques. In solving a problem, it finds a good compromise between the time efficiency typical of analytical solutions and the accuracy of the numerical methods, without the cumbersome analytical integration. Numerical results are presented for examples of catalytic and noncatalytic reactions, and the results are compared with the analytical results. The polynomial approximation technique provides results which are in good agreement with the numerical results and are more accurate than the flame sheet approximation results. The configuration of the system consists of a porous medium of total length Lt. Only the axial dimension x is considered. A premixed gas stream is fed into the system at x = -LIand flows toward the reaction front that is located at the position x = 0 inside the porous medium. Three zones are identified within the porous medium. Zone I is the preheating section. The reaction occurs only in zone 11, and the postreaction section is called zone 111.
0 1995 American Chemical Society
Ind. Eng. Chem. Res., Vol. 34,No. 3, 1995 795
2. Overview of the Methods 2.1. Flame Sheet Approximation (FS).When the activation energy is large, the reaction zone will be thin and at the scale of the bed length its thickness could be neglected (Matkowsky and Sivashinsky, 1979). To analyze the flame zone, the dependent variables of the system (temperature, concentration of species, etc.) are expanded in the form of perturbation series; i.e., for the temperature,
1
T,=T,+-T,+
Y
...
E R Tfo
Furthermore, we will use y to stretch the axial distance in zone I1 as follows: u
(3)
=xy
Although the reaction zone thickness is neglected at the actual scale x, it is assumed to be relatively large a t the scale u. The series expansions are substituted in the differential equations that describe the model. Then, by matching terms of the same order in y it will produce a larger set of new equations which hopefully will be easier to solve than the original ones (if y is large only the first perturbation variable will need to be evaluated). The integrated equations will relate the "limiting" values of the dependent variables of the adjacent nonreaction zones. A detailed account of the steps required is given elsewhere (Escobedo and Viljoen, 1993). 2.2. Polynomial Approximation Method (PA). This technique is proposed as a semianalytical approach to integrate a single ordinary differential equation or a set of coupled equations. It is based on the concept of spline functions and the discretization of the integration path. Suppose an unknown function y = y(x), defined by an nth order ordinary differential equation (ODE),is to be integrated by dividing the integration interval (XOJk) into 12 subintervals. Assume that the differential equation is satisfied at each node including the boundaries:
zdx"-) (k + d"-'yi d"yi
,...,
= 0,
1 equations) (4)
The boundary conditions at the ends of the interval are written as
2)
go(xoJo,% ,...,
gk Xk,yk,x, = 0, (n - n, equations) dyk dn-lyk) .e.,=
+
xifl = xi hi, (k equations)
(8)
L, = xk - xo, ( 1 equation)
(9)
+ +
The total number of equations is 312 n 1. The total number of unknowns, including yi, dyildx, ..., d"yildxn ( n 1for each of k 1points) and the Axi and xi (k of each) is ( n l)(k 1) 212. The degrees of freedom found by subtracting the number of equations from the number of unknowns is nk. This means that, in order to solve the system, we need t o know n additional equations for each of the k steps. There are several ways to provide these additional equations. Some of these ways will be discussed next, keeping in mind that different ways can be applied at different subintervals in the same problem: (i) If an analytical solution or even an approximate analytical solution is available for a given subinterval, then
+
+ + +
fi(%y,C,,C,,...,C,) = 0
(6)
The positions of the nodes (or step sizes) are specified by arbitrary equations like
(10)
where the Cis are the n integration constants (undetermined at this point). Equation 10 can be used to generate y, dyldx, ...,d"y1dx" by successive differentiation and applying the resulting expressions to the subinterval points i and i 1 to produce 2(n 1)equations. However, since the ODE (eq 4)is satisfied by eq 10, only 2(n 1) - 2 = 2n equations are independent. Subtracting the n unknown integration constants, we are left with n net additional equations as needed. (ii) A polynomial P(x) of order m (unknown at this point) is assumed to describe the function y(x) in s consecutive subintervals. At the first point i, it is assumed that up to the order rl the derivative of P ( x )is equal t o the corresponding derivative of the ODE (i.e., yi = P ( x ~dyildx ), = d P ( ~ i ) / d x..., , d'lyildx'l = d'lP(xi)/d3cr1). The corresponding top order derivatives that satisfy analogous conditions for the s - 1 interior points and the end point are rz and 7-3, respectively. Also assume also that r1, 7-2, and 7-3 I n. Although derivatives of higher order than n could be necessary at a given point (i.e., if rl > n at point i), no new variables will be considered. For each one of them an independent defining equation can be generated by repeated differentiation of the ODE (eq 4). It is necessary that
+
+
+
no. of eqs = s x n
+ no. of unknown coeff of polynomial
Hence,
(rl
+ 1 ) + (r3+ 1 ) + (s - l ) ( r z+ 1 ) = sn + m + 1 (11)
= 0, n, (n, equations, n, < n )
(5)
i
(k - 1 equations) (7)
We also have
+
and the inverse of the dimensionless activation energy (7) is used as the perturbation parameter. The definition of y is
Y=-
= 0,
This relationship allows us to determine m. The most important case to be considered here is when r1 = rz = 7-3 = n; in that case from eq 11,
m=n+s
(12)
The polynomials corresponding to this case will be identified by a double subindex as PA(n,s). For ex-
796 Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995
ample, for s = 1 we can derive the following relationships:
and
1 &Yi c , = -J
j!
ax’’ j = 0, 1,...,n
-+ n
1
(14)
For the parabola PA(1,l) the n = 1additional equation per interval required follows from eqs 13 and 14: (15) which is a very well-known approximate integration formula. The two additional equations for the cubic PA(2,l) are
Equations analogous t o eqs 13 and 14 can be derived for s = 2. The additional equation for the cubic PA(1,2) fitted to the points (i = 1, i and i 1)is
+
60
where
r=
+r)
dx
(%r,
q = f l ’ j = $1 q
- 1)
(19)
The two additional equations for the fourth order polynomial PA(2,2) are
d2Yi 2 d2Yi-q (1 r)(3 r ) -- r -
+
+
dx2
dx2
- 1 interior points satisfy the ODE (eq 4) while the outmost points satisfy continuity conditions of the function and first derivative. The analysis of this case can be found elsewhere (Villadsen and Michelsen, 1978). It only needs to be pointed out that it is possible to combine any piecewise approximate polynomial solution with an analytical one and maintain regularity of the composite solution if care is taken in accounting for all the degrees of freedom of the system as shown before. 2.2.1. Discussion. Expressions like eqs 15-21 could also be generated by suitable truncated Taylor expansions. These equations are implicit integration formulas and are related to the predictor-corrector integration approach, although no predictor formula is relevant for our purposes. In general, the nodes can be determined by any adequate procedure including orthogonal collocation. It is beyond the scope of this work to explore the capabilities of the PA(n,s) as compared to other conventional numerical integration techniques. Our interest lies in the lowest polynomial orders and the fewest number of subintervals. If in a given application the number of steps which require numerical integration can be kept very low, i.e., less than 2, then the moderate number of resulting algebraic equations t o be solved would make the approach comparable t o an analytical technique. The computational time required would be substantially less than for a truly numerical integration, especially in boundary value problems. The piecewise combination with an analytical solution can produce some computational advantages by resolving the stiffness of the ODE. In this context the PA(n,s) discussed earlier and in particular the PA(n,l) and PA(n,w) are advantageous since they take full advantage of the information available at the step points. When approximate analytical solutions are used in a given subinterval (case i), the differential equations may not be exactly satisfied at either end of the subinterval. Therefore no continuity of derivatives is assured (loss of regularity). Also note that, in general, the step sizes are unknown a priori and in some cases they can be relatively large. In those cases, the PA(n,l) is especially advantageous since for solving one step information only the initial and final point is used. Other integration formulas available (of the predictor-corrector type) require more information from previous points. Explicit integration techniques are less stable (Le., direct Taylor series expansions) and since the step size is unknown a priori, methods such as the fourth order Runge-Kutta scheme would be more complex t o apply.
3. The Models. General Assumptions
(21) where r, q , and j are again given by eq 19. (iii) A situation exists similar to case ii but now rl, r2, and r3 can be less than n. In this case some of the original variables cPyf&j and some of the eqs 4 can become meaningless (not used) and the analysis of degrees of freedom would need to be reformulated. An especially interesting case arises when rl = 1-3 = 1and r2 = n since it would correspond to a typical orthogonal collocation in finite elements scheme where only the s
Descriptions of the models are presented next, followed by the alternative solutions. The following assumptions have been adopted for both the one-phase and two-phase models: 1. Radial and angular variations in temperature or species concentration are neglected. 2. The wall of the bed is adiabatic. 3. The reaction is highly exothermic with a relatively large activation energy. 4. For simplicity’s sake, a first order global reaction rate will be used. Other typical kinetics can be handled similarly. 5. The change in molar flux along the bed is neglected (equimolar reactions or large excess of inert).
Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995 797
6. Axial mass diffusion effects are neglected because mass convection is the dominant mode of transfer (Gatica et al., 1987). 7. The pressure drop over the bed is neglected because the bed is shallow (i.e., less than 0.20 m) and the gas velocity is usually low (Sathe et al., 1990; Gatica et al., 1987). 8. The ideal gas law is used as the equation of state. 9. All physical properties are considered constant along the bed. 10. Net heat losses through the inlet cross section are neglected. 11. It is assumed that the whole outlet cross section can radiate heat to the load. This means that the net radiation of the solid surfaces close to the outlet and in the direct line of sight of the outer space is equivalent t o the radiation of a surface of area equal to the bed cross section and temperature equal t o the solid phase at the outlet.
4.1. Analysis of Nonreacting Zones. The solution of eq 22 when the reaction term is dropped can be expressed as
where c = GC,
(30)
T j denotes the temperature at the boundary with the reaction zone. For zone I
T j = Tl,
z = -Tin
(31)
The slope at the interface of zones I and I1 would be
4. One-Phase Model
When the heat transfer coefficient and transfer area between the fluid phase and the solid matrix are large, thermal equilibrium can be assumed between both phases which corresponds to the pseudo-homogeneous or one-phase model (Froment and Bischoff, 1979). Energy and component balances yield d
-[(k, dx
+ b P ) $d
dT
+
- GCp dx ( - W R =0 dw-
R
dx
G
(22)
k,=-
-EIRT
(24)
ckJI
R
+ b P ) d T = GCp(T - Tin) w = win
4.2. Solution of Reaction Zone by FS Approach. A slight variation of the analysis presented by Escobedo and Viljoen (1993) is summarized next. When the expanded variables (cf. eqs 1-3 are substituted in eqs 22-24) and terms of the same order are considered, it is found that Tm = T, (the maximum temperature) and the matching conditions for the reaction zone are at v -m , (T = TI, w = Win)
4
The radiation effect is described by the diffusion approximation (cf. Norbury and Stuart (1989) and Vortmeyer (1980)). The use of the diffusion model to describe radiation is acceptable as long as the medium can be considered as optically dense (Sparrow and Cess, 1970) as it is usually found in catalytic beds ( E > 0.7). Otherwise, alternative models should be used (Viskanta and Menguc, 1989). Since the use of the alternative models would make a n analytical treatment extremely complex and much less appealing, the diffusion approximation will be retained at the cost of sacrificing some accuracy in the case of porous burners. On the other hand, the diffusion model of radiation can produce approximate results for systems with a porosity as large as 0.947 (Vortmeyer, 1980). The boundary conditions are, at the entrance, x = -L1,
(k,
The slope at the interface of zones I1 and I11 would be
(23)
where for a first order kinetics and ideal gas behavior the reaction rate is given by
R = k,-ew T
For zone 11,
(26)
(27)
T, = y(Tl - T,) when v
-
+m,
and
dT,
(T = Tk, w = wout)
The differential equations of O(y) are dTl &
dTk & - P i n - Bout
where Pin,out
- (-W f l - k, + bTm3uwin,out
and a t the outlet (k,
+bp)
= h,(T$ - T?
(28)
a=
ck,IIe-EIRTm
GE
dT,
=-
dx
798 Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995
After integrating eq 38 between u and -=, we obtain
eq 22 and the result integrated from x = -L1 t o x , we get, at any point,
IBy taking the limit u
-
+w
Equations 45 and 46 can be combined to produce an expression of the form
in eq 41, we get
(47)
At u = 0, T = Tm and dTnIdv = 0; eq 41 yields
At this point TI and Tk are undetermined. As the thickness of the reaction zone is assumed to be small when the flame sheet approximation is used, the differences between TI, Tk, and T, tend to vanish. Therefore, we can neglect these differences when computing the outer slopes from eqs 32 and 34. However these temperature differences cannot be neglected in the arguments of the exponentials which appear in eqs 42 and 43. If y is a very large number (due to a large activation energy), the exponential terms tend to zero even for a small temperature difference. These exponential terms can be finite if the reaction zone actually extends to one end of the bed because T, might be marginally larger than the temperature a t the end. Therefore, TOand Tk+l will be used instead of TI and Tk, respectively, in the exponential terms of eqs 42 and 43. When these approximations are adopted and dTlldx from eq 33 and the definitions of a and Bin are substituted in eq 43, the following expression can be obtained:
r
Therefore, a cubic polynomial will be fitted to each subinterval. The analytical solutions are now available for the first and last integration steps. However, for all practical purposes these analytical solutions play the same role as the approximate solutions given by the polynomials of the PA method (i.e., the differential equations are satisfied at the step end points). Therefore, in order to solve the system, the step points 1 and k have to be specified (cf. eq 7). There is no obvious way of specifying the limits of the reaction zone, because the reaction rate never becomes identically zero. But we can borrow the concept of terms of similar order from our perturbation series method and specify that the reaction rate at the ends of the reaction zone is 1order of magnitude smaller than at the middle point, i.e., dw,
-=u
dw, ~
dx
dwkand --
dx
dw2
UkdX
(48)
-
where u1 and U k are chosen to be O(1ly) 0.1. Equations 16,17,47, and 48, (fori = 1,2 and having A x 2 = A x g ) along with the boundary conditions of the nonreaction zones, eqs 32 and 34, close the system. This system of equations is easy to solve without resorting to complex routines for simultaneous nonlinear equations. Equation 48 is meaningless if either zone I or I1 or both do not exist. This inconsistency will become apparent from the calculated results. The way to correct it is to replace the faulty equatiods) by the corresponding boundary condition(s).
5. Two-Phase Model
Equations 32, 34, 37, 42, and 44 relate the unknowns Tm, dTlIdx, dTkldx, Tk+l, TO,and wout. L2 is evaluated from eq 29 applied a t T - Tk+l (zone III), and TOis evaluated from eq 29 applied at x = -L1= Lt - La (zone I). This system of equations can be easily solved iteratively. 4.3. Solution of Reaction Zone by PA. To set some ideas, the reaction zone is divided into two sectors (points 1-3) and point 2 is halfway between the beginning and the end of the reaction zone. Applying the differential equations of the model to each point i, we have
['
d2Ti=
h2 kei
JL) + c dT, - (-AE?-e k*wi
-3b( Ti dTi 2
Ti
-E/RT,
1
(45)
where k,i = k, + bTi3. If R from eq 23 is substituted in
5.1. Model Description. In this case no thermal equilibrium exists between the gas and solid phase. The additional assumptions for this model are as follows: 1. The internal radiation effects will be neglected. Although this assumption can be a significant limitation if the system reaches very high temperatures as in porous burners, a qualitative correction can be made by speciwng a larger effective solid conductivity. This assumption is necessary to keep the analysis tractable. 2. There is negligible gas phase conduction (Singh et al., 1991; Viskanta, 1991). 3. For a porous burner the reaction will be assumed to take place strictly in the gas phase (the porous matrix is noncatalytic). 4. For a catalytic bed the reaction will be assumed to take place inside the solid pellets and therefore the heat is released in the solid phase. In this case a global reaction rate model is assumed which may take into account the film resistance due to mass transfer from the bulk gas stream t o the pellet surface and internal pore diffusion efficiency factors.
Ind. Eng. Chem. Res., Vol. 34,No. 3, 1995 799 differential energy balance for the solid:
d2T, ke- h&T,- Tg)- Y(-AWR,
dx2
(N - M), (49)
differential energy balance for the gas:
At x = 0 (interface of zones I and 11) we have
+
dTg h,(Ts - Tg)+ (1- Y)(-AH)Rg = 0 -GC pax
Ts1= Tin + fA + fnB
(65)
dTgl = \Tsl d x c
Tgl)
(66)
- q T g l - Tin) ke
(67)
(50) The differential component balance and reaction terms are given by
-dTsl
dx R, = k,-ew T,
-EJRTS
and R~ = k,oe-E/RTg
Tg
(52)
+
ho Tg= Tin -(T, GCP
+ T, = Ti, + (Tgl- Tin)fmeMx Tg= Tin (Tgl- Tin)eMx
- Tin)
C=T
(54)
(55)
dx
dT, ke= h , ( c - T') h,(Tg- T,)
+
+ AeMx+ BeNx T, = C + AfmeMx+ BfneNx
+%7,
M >N
f,=l+-, cM f n = l + cN-
h,
h,
dTgk
dTsk
dx
dx
-c-+k$Mc(M - N)
(72) At the interface of zones I1 and I11 we have Tg= Tgk,T, = Tsk,and (73)
(57) (58)
The boundary condition at the outlet is written as
k,--
dTs(k+lj
dx
- hc(Tg(k+l) - T s ( k + l ) ) - hr(Ts(k+l:
-
5":) (74)
(59) where (60)
Remark: The numerical values of M and N will be the same for both zones if the physical properties involved are unchanged. For zone I
C = Tin
dx ; B =
(56)
5.2. Analysis of the Nonreaction Zones. For zones I and I1 the solution has the common form
Tg= C
c(M - N)
(71)
c dx
dTsk
dTgk
A=
* (1
(70)
ke dTsk gk
and at the outlet
M , N = &2c [-i
(69)
For zone I11 we have
c--ka-
where
(68)
T,, = Tin+ fm(Tgl- Tin)
(53)
w = Win
dx
-
If h, >> c as is usually the case, then N
