Ind. Eng. Chem. Res. 1999, 38, 767-777
767
Analytical Criteria for Validity of Pseudohomogeneous Models of Packed-Bed Catalytic Reactors Sandra M. S. Dommeti and Vemuri Balakotaiah* Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4792
David H. West The Dow Chemical Company, Freeport, Texas 77541
Pseudohomogeneous models are used extensively to simulate the behavior of packed-bed catalytic reactors. The use of highly active catalysts increases the importance of the transport resistances between the solid and the fluid, and hence, may invalidate the applicability of the pseudohomogeneous models. There exist criteria in the literature that predict when pseudohomogeneous models can be used to approximate the behavior of catalytic reactors. These criteria are derived by assuming that the predictions of pseudohomogeneous and the two-phase models differ only by a small amount and hence Taylor expansions can be used to quantify the difference between the two models. In this work, we show that the existing literature criteria may lead to erroneous conclusions and present new criteria that may be used to determine the regions of validity of pseudohomogeneous models. For example, when Lep ) h/FfCpfkc g 1, we show that, for a single exothermic first-order reaction, the pseudohomogeneous model has the same qualitative features as the two-phase model, if k(T0)/(kcav) < exp{4Lep - 2 - [E/(RT0)][(-∆H)c0]/[FfCpfT0]}, where Lep is the particle Lewis number, Ff is the density of the fluid phase, Cpf is the specific heat at constant pressure of the fluid phase, k(T0) is the reaction rate constant at the feed temperature T0, R is the universal gas constant, E is the activation energy, av is the external particle surface area per unit volume of the reactor, kc is the mass-transfer coefficient, h is the fluid-solid heattransfer coefficient, (∆H) is the heat of the reaction, and c0 is the feed concentration. A similar criterion is presented for 0 < Lep < 1. These criteria are very different from those in the literature and are derived from a rigorous analysis of the single-phase and two-phase models and comparison of their qualitative features (bifurcation diagrams). Application of the criterion to typical industrial situations show that they could differ from that of Mears (Mears, D. E. J. Catal. 1971, 20, 127),1 by a factor of 10-104. It is also shown that the two-phase plug flow model has 10 qualitatively different types of bifurcation diagrams. Some of these contain hightemperature branches on which the solid temperature exceeds the adiabatic temperature rise. 1. Introduction Mathematical models describing the steady-state and dynamic behavior of packed-bed catalytic reactors are obtained by writing down the continuity, momentum, species, and energy balances and combining these with the various constitutive equations for inter- and intraphase heat- and mass-transfer and chemical reaction processes. Depending on the level of detail included, these models can vary from a pair of ordinary differential equations to several partial differential equations in two/three spatial coordinates and time. In addition, the model equations describing catalytic reactors are highly nonlinear because of the strong coupling between the transport processes and the chemical kinetics and the nonlinear dependence of the kinetic and transport coefficients on the state variables (temperature or concentrations). The model equations also contain several parameters, which are known or can be estimated with uncertainties ranging from 0 to 20% or 30%. It is well-known that these nonlinear models exhibit multiple steady states, periodic states, and various spatiotemporal patterns. * Corresponding author.
From a practical point of view (for design and control purposes), it is desirable to have “simpler models” that retain all the qualitative features of the system. However, how to determine such models is not at all obvious without analyzing hierarchies of models. In the literature, the pseudohomogeneous models are developed by assuming that the interphase transport processes are much faster compared to the reaction processes. Criteria for the validity of these models are developed by assuming that these model predictions differ from the two-phase model predictions only by a small amount and hence Taylor expansions can be used to quantify this difference. While two-phase models of catalytic reactors were developed by Wicke,2 Liu and Amundson,3 Luss and Amundson,4 and many others, they were not analyzed in any detail. This was mainly due to the lack of mathematical techniques such as singularity and bifurcation theories (which were developed in the 1970s and 80s) that are extremely helpful in determining the qualitative features of the models and in guiding the computer simulations. The analysis of the simplest twophase model, that is, a two-phase CSTR, has been done only recently5 and has shown several surprising new results such as the existence of isolated high-tempera-
10.1021/ie980365g CCC: $18.00 © 1999 American Chemical Society Published on Web 01/20/1999
768
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
ture branches (with solid temperatures exceeding the adiabatic temperature) and oscillatory behavior in adiabatic reactors. The main purpose of this article is to examine the qualitative features of two-phase catalytic reactor models and develop criteria under which they simplify to the so-called pseudohomogeneous models. We examine the simplest two-phase model of a tubular catalytic reactor that accounts for the interphase heat- and masstransfer resistances, that is, the plug flow model with interphase resistances. We determine the different possible types of steady-state bifurcation diagrams (of solid-phase temperature versus the Damko¨hler number) and classify the steady-state behavior in the parameter space. These results are then used to determine the conditions under which the two-phase model behavior is the same as that of the pseudohomogeneous model (which neglects the interphase heat- and mass-transfer resistances). In the next section, we formulate a general onedimensional two-phase model of a catalytic reactor and review some limiting cases of this model. In section 3, we review briefly the recent results of Subramanian and Balakotaiah5 on the steady-state behavior of the twophase well-mixed reactor model. In section 4, we analyze the two-phase plug flow model and present some new results. In section 5, we compare the qualitative features of the two-phase models with the pseudohomogeneous models and arrive at analytical criteria for the validity of the latter. In the last two sections, we compare the new criteria to the literature criteria, interpret the practical relevance of the results of this work, and discuss some possible extensions.
∂θs ) 0, ∂ξ
@ξ ) 0, 1
∂xf ∂θf ) 0, ) 0, ∂ξ ∂ξ
(5)
@ξ ) 0, 1
(6)
xf 1 ∂xf ) 0, 2 ∂ξ Da φm
@ξ ) 0
(7)
θf ∂θf ) 0, 2 ∂ξ Da φh
@ξ ) 0
(8)
where
t ) k(T0)t′, ξ )
(
)
Ti - T 0 z , θi ) γ L’ T0 c0 - ci , xi ) c0
i ) s, f
λs k(T0) k(T0)FfCpf , Dapm ) , Daph ) , λf kc a v hav Dapm FsCps E h , γ) , Lep ) ) , σ) FfCpf RT0 Daph FfCpfkc MmCpm (-∆H)c0 Le ) 1 + , , B)γ FfCpfT0 VR[FfCpf + (1 - )FsCps]
Λ)
Da )
L2k(T0)FfCpf L2k(T0) k(T0)L , φh2 ) , φm2 ) u λf De
and 2. Mathematical Models of Catalytic Reactors A mathematical model for an adiabatic tubular catalytic reactor has been developed by Balakotaiah et al.6 The model is based on the assumptions of negligible temperature and concentration gradients inside the particles, flat velocity profile, negligible radial gradients, and constant catalyst activity. For the case of a firstorder exothermic reaction, the model equations in dimensionless form may be written as 2 ∂xf 1 ∂ xf 1 ∂xf xs - xf ) 2 2+ ∂t Da ∂ξ Dapm φ ∂ξ m
(1 - )
L
(1)
( )
xs - xf ∂xs )+ (1 - xs) exp ∂t Dapm
θs
θs 1+ γ
2 ∂θf ∂ θf 1 ∂θf θs - θf ) 2 2 + ∂t Da ∂ξ Daph φ ∂ξ
(2)
(3)
h
and
∂θs (1 - )Λ ∂2θs θs - θf ) + σ(1 - ) ∂t Daph φ 2 ∂ξ2 h
( )
B(1 - xs) exp
with boundary conditions,
θs
θs 1+ γ
(4)
L ) Le + (Le - 1)σ(1 - )
(9)
Here, xf, xs, θf, and θs represent the dimensionless conversions and temperatures in the fluid and solid phases, respectively. The above model is written in terms of the mass and thermal Thiele moduli (φm2 and φh2) instead of the mass and heat Peclet numbers (Pem ) φm2/Da, Peh ) φh2/Da) in order to isolate the residence time in the Damko¨hler number. Another advantage of the above formalism is that various limiting cases of the model can be seen more easily. The parameters of the model are the porosity (), the ratio of the solid-to-fluid thermal conductivities (Λ) and heat capacities (σ), the reactor Lewis number (Le), the reactor Thiele modulus (φm2) which is the ratio of the longitudinal dispersion time to the reaction time, the thermal Thiele modulus (φh2), which is the ratio of the longitudinal thermal dispersion time to the reaction time, the Damko¨hler number (Da), which is the ratio of the convection time to the reaction time, the dimensionless activation energy (γ), the dimensionless adiabatic temperature rise (B), the particle mass Damko¨hler number (Dapm), which is the ratio of the interphase mass-transfer time to the reaction time, and the particle heat Damko¨hler number (Daph), which is the ratio of the interphase heat-transfer time to the reaction time. The ratio of Dapm to Daph is the particle Lewis number (Lep). We note that the model defined by eqs 1-8 reduces to the more commonly used pseudohomogeneous model in the limit of vanishingly small particle heat and mass Damko¨hler numbers (Daph f 0 and Dapm f 0 with all other parameters fixed). In this limit, the concentration
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 769
and temperature in the two phases are nearly identical at all points along the reactor and eqs 1-4 may be combined to give the pseudohomogeneous model defined by the following equations:
The steady-state well-mixed two-phase model is obtained by dropping the time derivative terms in eqs 1518 and is given by the algebraic equations
-
∂θ Le* ) ∂t
( )
1 ∂θ 1 ∂2 θ θ + B(1 - x) exp (10) 2 2 Da ∂ξ θ φh* ∂ξ 1+ γ
-
( )
m
θ
1+
θ γ
(11)
θs
θs 1+ γ
-
Le* ) Le + σ(1 - )Le and FfCpfL k(T0) [λf + (1 - )λs]
(12)
1 ∂x x ) 0, @ξ ) 0 2 ∂ξ Da φm
(13)
and
θf ) Bxf
(14)
The above pseudohomogeneous model may be reduced further to the homogeneous model in the limit λs ) λf (φh*2 ) φh2) and FfCpf ) FsCps (σ ) 1 or Le* ) Le). It also reduces to the homogeneous model in the more obvious limit of f 1. Various other limiting cases of the model defined by eqs 1-8 may be formulated but we consider here only two special cases. The first limiting case is that of a wellmixed two-phase system obtained in the limit φh2 f 0, φm2 f 0 (with Λ finite). In this case, the model simplifies to
xf xs - xf dxf )+ dt Da Dapm
(15)
( )
dxs x s - xf (1 - ) )+ (1 - xs) exp dt Dapm
L and
σ(1 - )
θs
1+
θs γ
θf θs - θf dθf )+ dt Da Daph
xs )
θs -
(1 - )
Lep (Dapm + Da) θ B (Dapm + LepDa) s
(25)
( )
θs 1+ γ
]
( ) θs
θs 1+ γ
) 0 (26)
∂xf 1 ∂xf xs - xf )+ ∂t Da ∂ξ Dapm
(27)
( )
xs - xf ∂xs )+ (1 - xs) exp ∂t Dapm
(16) L
(18)
(23)
( )
[
σ(1 - )
θs 1+ γ
(22)
The second limiting case of interest is the two-phase plug flow model obtained in the limit φh2 f ∞ and φm2 f ∞. In this case, the model equations simplify to
and
θs
)0
θs
B (Dapm + LepDa) 1 Lep Lep (Dapm + Da) θ exp B (Dapm + LepDa) s
(17)
θs - θf dθs )+ B(1 - xs) exp dt Daph
θs 1+ γ
(24)
xf ) Da(1 - xs) exp
and
∂x ∂θ ) 0, ) 0, @ξ ) 1 ∂ξ ∂ξ
θs
Combining eqs 19-22, we get
The boundary conditions are
1 ∂θ θ ) 0, @ξ ) 0 2 ∂ξ Da φh*
(20)
(21)
( )
θs - θf + B(1 - xs) exp Daph
2
φh*2 )
)0
θs - θf θf )0 + Da Daph
and
where
(19)
( )
x s - xf + (1 - xs) exp Dapm
and
1 ∂2 x ∂x 1 ∂x ) 2 2+ (1 - x) exp ∂t φ ∂ξ Da ∂ξ
xf xs - xf )0 + Da Dapm
θs
θs 1+ γ
∂θf 1 ∂θf θs - θf )+ ∂t Da ∂ξ Daph
(29)
( )
∂θs θs - θf )+ B(1 - xs) exp ∂t Daph
(28)
θs
θs 1+ γ
(30)
The steady-state version of this model is a differentialalgebraic system in terms of the solid and fluid temperatures:
770
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
dθf Lep (θ - θf), θf ) 0 @ Da ) 0 ) dDa Dapm s
( ) ( )] θs
Dapm(B - θf) exp
θs 1+ γ θs - θf ) θs Lep 1 + Dapm exp θs 1+ γ
[
(31)
(32)
( ) ( )]
xs )
[
1 + Dapm exp
θs γ
(33)
θs
1+
θs γ
3. Review of Results for a Two-Phase Well-Mixed Model We review in this section the results obtained by Subramanian and Balakotaiah5 on the steady-state behavior of the two-phase well-mixed model. These results are needed for the understanding and interpretation of the behavior of the two-phase plug flow model presented in the next section. As shown in section 2, the steady-state behavior of the two-phase well-mixed model is given by a single algebraic equation for θs eq 26, which may be rewritten as
( )
θs exp Da )
θs
(
B - Dapm - θs θs Lep 1+ γ B - θs
) (34)
Subramanian and Balakotaiah5 applied singularity theory to analyze the different possible bifurcation diagrams of θs versus Da and showed that there are 10 such qualitatively different diagrams for the case of γ f ∞ (Figure 1). These diagrams are obtained for (B, Lep, Dap) values in the 10 regions of the phase diagram obtained by plotting the hysteresis, isola, and boundary limit varieties. The hysteresis variety is the locus of the parameter values where two limit points (an ignition and extinction point) come together and disappear (or appear). It is given by the expression
B)
4Lep 2
Lep + e Dapm(1 - Lep)
]
(36)
and exists only if Lep < 1. The boundary limit set is the locus of parameter values at which a limit point appears or disappears at Da ) 0. This locus may be expressed in a parametric form
Lepθs2 θs - 1
Dapm ) (θs - 1) exp(-θs),
1 < θs < ∞
(37)
When interphase gradients are ignored and a pseudohomogeneous model (Dapm ) Daph ) 0) is used, that is, eq 34 is simplified to
θs
1+
Lep
Dapm(1 - Lep)
B)
The corresponding conversion profiles are obtained by using eq 23 and the relation
xf + Dapm exp
[
B ) ln
θse-θs Da ) B - θs
(38)
only two types of bifurcation diagrams (single-valued and S-shaped) exist. Thus, the introduction of interphase gradients (which adds two extra parameters, namely Dapm and Lep) gives eight additional types of behavior. The most important differences between the two models are (i) the existence of high-temperature branches (for Lep < 1) in which the solid phase can attain temperatures as high as B/Lep which is higher than the adiabatic temperature rise B (the fluid temperature never exceeds this value) and (ii) the possibility of the highest solid-phase temperature being attained at zero residence time (Da ) 0). Subramanian and Balakotaiah5 have also analyzed the stability of the various solution branches in Figure 1 and showed that when the reacting fluid is a gas, the stability changes only at the limit points. Therefore, the middle solution branches are unstable while the upper and lower branches are stable (as shown in Figure 1). 4. Analysis of the Two-Phase Plug Flow Model For the case of γ f ∞ (positive exponential approximation), the two-phase plug flow model is described by the following differential-algebraic system:
(35)
The isola variety defines the locus of parameter values at which two branches of the bifurcation diagram intersect and separate. (This is called the hyperbolic isola variety in contrast to the elliptic isola variety at which a new solution branch emerges as an ellipse. The defining conditions for both are the same.) The isola locus of eq 34 may be shown to be
Figure 1. Schematic bifurcation diagrams of solid temperature versus Damko¨hler number for two-phase lumped and plug flow models.
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 771
dθf Lep (θ - θf), θf ) 0 @ Da ) 0 ) dDa Dapm s
(39)
(B - θf) exp(θs) θs - θf ) Dapm Lep[1 + Dapm exp(θs)]
(40)
When Dapm ) 0 (no interphase gradients), θs ) θf ) θ and we get the usual pseudohomogeneous plug flow model:
dθ ) (B - θ) exp(θ), θ ) 0 @ Da ) 0 dDa
(41)
for which there is always a unique solution (i.e., bifurcation diagram of θ versus Da is single-valued and monotonic). When Dapm > 0, the derivative of eq 40 w.r.t. θs can vanish and hence eqs 39 and 40 define a differential-algebraic system with an index greater than unity.8 To our knowledge, the bifurcation behavior of such systems has not been studied in the literature. Since the particle eq 40 can have multiple solutions and the plug flow model allows discontinuous solutions (due to the absence of conduction in either the solid or fluid phase), the two-phase plug flow model defined by eqs 39 and 40 can have infinite number of solutions. This was pointed out by Wicke,2 Liu and Amundson,3 and others. Eigenberger9 showed that when conduction is included in the solid-phase energy balance, the number of solutions is finite. Balakotaiah10 has shown that many of these solutions emerge from the unstable middle branch are nonphysical (with sinusoidal temperature profiles) and hence are unstable. We attempt to show here that when only physically meaningful solutions are considered, the behavior of the above model is qualitatively similar to that of the two-phase well-mixed system reviewed in the previous section. We note that when θf ) 0 and Dapmeθs . 1, eq 40 simplifies to
θs ≈
B Lep
Thus, as in the case of the well-mixed system, the highest solid-phase temperature can be B/Lep (when Lep < 1) and is attained at zero residence time (or at the inlet of the reactor). Similarly, for Da f ∞, we have from eqs 39 and 40
Thus, both the solid and fluid temperatures reach the adiabatic temperature rise when the residence time is very large (or at the exit of the reactor). It follows from eqs 39 and 40 that the fluid-phase temperature (θf) is a monotonically increasing function of Da and cannot exceed B. In contrast, the solid-phase temperature (θs) can be either increasing or decreasing with Da (or along the reactor length). We note that when Da ) 0 (or θf ) 0), the particle heat-balance equation given by (40) simplifies to
BDapm exp(θs) Lep[1 + Dapm exp(θs)]
f)
∂f ∂2f ) )0 ∂θs ∂θ 2
(42)
The solution of eq 42 determines the solid-phase temperature at the inlet. This equation can have three solutions for some range of Dapm values when B > 4Lep.
(43)
s
where
f ) θs - θf -
Dapm (B - θf) exp(θs) Lep [1 + Dapm exp(θs)]
Solving eq 43 gives the hysteresis locus to be
( )
B ) 4Lep - 2 + ln
1 Dapm
(44)
θf ) B - 4Lep
(45)
θs ) θf + 2
(46)
The Da value at the hysteresis point is obtained by integrating eqs 39 and 40 with the constraints given by eqs 45 and 46. In addition to the above loci, a bifurcation diagram of θs versus Da may change qualitatively when the parameters cross the isola locus defined by the equations
f)
∂f )0 ∂θs
df ∂f dθf ∂f dθs ) + )0 dDa ∂θf dDa ∂θs dDa
(47)
(48)
Combining eqs 47 and 48, we see that the last condition is equivalent to the equation
∂f )0 ∂θf
θs ) θf ) B
θs )
We define the boundary limit set as the locus of B, Dapm, and Lep values at which the number of solutions of eq 42 changes from one to three. This locus is given by eq 37. We note that when the particle equation (40) has a unique solution for θs for all values of θf g 0, B, Dapm, and Lep, then the two-phase plug flow model also has a unique solution. The transition from unique to multiple solutions occurs when the particle equation (40) begins to have multiple solutions. We define this locus to be the hysteresis variety of the two-phase plug flow model. It satisfies the equations
(49)
Elimination of θs and θf from eqs 47 and 49 gives the isola locus defined by eq 36. As in the case of the well-mixed two-phase model there are three different cross sections of the phase diagrams in the (B, Dapm) space. These will be discussed in more detail now. When Lep g 1, only the hysteresis and boundary limit varieties exist and divide the (B, Dapm) plane into three regions as shown in Figure 2 for a typical case of Lep ) 2. For (B, Dapm) values in region (i) the particle equation has always a unique solution and hence the bifurcation diagram (θs versus Da or θf versus Da) of the two-phase plug flow model is singlevalued and monotonically increasing. When (B, Dapm) values cross the hysteresis locus, the particle equation has multiple solutions for some range of θf (and hence Da) values. Here, the low-temperature branch that started at Da ) 0 ceases to exist at some critical Da
772
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
Figure 2. Steady-state classification diagram of the two-phase plug flow model for Lep ) 2.0.
(ignition point). Similarly, the high-temperature branch (for which θs ) θf ≈ B for Da f ∞) ceases to exist for Da values below a critical value (extinction point). Thus, the bifurcation diagram in region (viii) is S-shaped. A numerically computed bifurcation diagram of this type is shown in Figure 3. We note that for Da values exceeding the extinction value, an infinite number of profiles are possible since the plug flow model allows the jump from a low- to high-temperature branch at any value of Da between the two limits. When the (B, Dapm) values cross the extinction branch of the boundary limit set into region (x), the particle equation has three solutions (θsl, θsm, and θsh) even at Da ) 0. The lowtemperature solution starting at θsl continues until some critical Da (ignition point). Similarly, the high-temperature solution branch starting at θsh continues to θs ) B at Da ) ∞. When we tried to compute a middle unstable branch starting at θsm, it was found that it has a limit point at a Da value smaller than the ignition point (i.e., it does not join the lower branch at the same Da value (Figure 3)). We believe that this feature is unique to the plug flow model and disappears when conduction is included either in the fluid- or solid-phase energy balance. This observation is supported by some numerical calculations performed by Nguyen and Balakotaiah7 with the conduction term included in the solidphase energy balance. We note that the solid and fluid temperatures at the limit points are uniquely determined by the particle eq 40 and its derivative w.r.t. θs. The Da value at the limit point is obtained by integrating the differential-algebraic system defined by eqs 39 and 40. When, the particle eq 40 has multiple solutions at Da ) 0, the middle solution branch which has a higher temperature gives a smaller Da value at the limit point. In contrast, the lower temperature branch gives a higher Da value. However, when conduction is included, both these values coincide since the temperature gradient for the middle branch is negative while it is positive for the lower branch. Thus, the discontinuity in the Da value is a limitation of the plug flow model. If only the stable branches are considered, the plug flow model predicts the same qualitative behavior as that of more complicated models that include the conduction term in the energy balance.
Figure 3. Numerically computed bifurcation diagrams of the twophase plug flow model for Lep ) 2.0 and (B, Dapm) values in regions (viii) (top) and (x) (bottom), of Figure 2.
Figure 4. Steady-state classification diagram of the two-phase plug flow model for 0.5 < Lep < 1. Top figure: Schematic phase diagram. Bottom figure: Numerically computed phase diagram for Lep ) 0.6.
Finally, we note that when the parameters (B, Dapm) cross the ignition branch of the boundary limit set, the particle equation has again a single solution for all θf g 0 and the bifurcation diagram is again single-valued and monotonic. The typical phase diagrams for Lep values in the ranges 1/2 < Lep < 1 and 0 < Lep < 1/2 are shown in Figures 4 and 5, respectively. In Figure 4, the hysteresis, isola, and boundary limit sets divide the (B, Dapm) space into seven regions, each corresponding to a dif-
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 773
Figure 5. Steady-state classification diagram of the two-phase plug flow model for Lep ) 0.1.
ferent type of bifurcation diagram. Since some of these regions are very small and are not visible in the numerically computed phase diagram for Lep ) 0.6, we have shown a schematic phase diagram at the top of Figure 4. As the transitions that occur when the parameters cross the hysteresis and boundary limit varieties are already described, we consider only the isola locus. For parameter values on this locus, θs ) B is a solution for all Da. This branch may intersect other solution branches emanating from Da ) 0. When the parameters cross the isola locus, the intersecting branches are separated, giving rise to solution curves either above the line θs ) B or below the line θs ) B or both. Figure 5 shows the phase diagram for Lep ) 0.1. Once again the varieties divide the parameter space into seven regions corresponding to seven different types of bifurcation diagrams. Four of these diagrams also exist for 1/2 < Lep < 1 and three are new. All of the 10 different diagrams are generated for 0 < Lep < 1. We show in Figures 6-8 the numerically computed bifurcation diagrams for some typical parameter values. It should be pointed out that the classification given above may be incomplete because of the existence of additional solution branches (stable or unstable) at intermediate values of Da. By analyzing eqs 39 and 40, we can prove that a bifurcation diagram of θs versus Da cannot have more than three limit points. When three limit points exist, the qualitative nature of the bifurcation diagram changes when the relative position of these limit points changes. The numerical calculation shown in the top diagram of Figure 7 indicates that a double-limit variety exists. The classification given here can be completed by plotting the double-limit variety, but this will not be pursued here. It is interesting to note that when Lep < 1/2, the plug flow model predicts three stable steady states (top diagram of Figure 7). One of these is the quenched steady state, and the second one is the ignited steady-state with the temperature close to the adiabatic temperature. The third and new stable steady state is the ignited state with θs between B and B/Lep.
Figure 6. Numerically computed bifurcation diagrams of twophase plug flow model for Lep ) 0.8 and (B, Dapm) values in regions (ii) (top) and (ix) (bottom).
Figure 7. Numerically computed bifurcation diagrams of the twophase plug flow model for parameter values in regions (iii) (top) and (iv) (bottom).
5. Analytical Criteria for Validity of a Pseudohomogeneous Model Having determined the qualitative features of the two-phase plug flow model, we can now examine the conditions under which the behavior of the heterogeneous model reduces to that of the pseudohomogeneous plug flow model. The latter has a unique solution and a single-valued bifurcation diagram of θ versus Da. It follows from Figure 2 that when Lep g 1, the two-phase model has a single-valued bifurcation diagram if the (B, Dapm) values are in region (i). Within the practical range of interest (0 < Dapm < e-2), the boundary of region (i)
774
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
the criterion given by eq 50 may be written as
Da/pm < e4Lep-2
(54)
We note that Da/pm increases exponentially with Lep. When 1/2 < Lep < 1, the boundary of region (i) is formed by the boundary limit and hysteresis varieties (Figure 4). For 0 < Lep < 1/2, this boundary consists of a segment of the extinction branch of the boundary limit set and a segment of the isola variety. We note that the isola points do not exist if
Dapm
