Singularity Theory Approach for Calculating the Runaway Boundaries

3230

Ind. Eng. Chem. Res. 1997, 36, 3230-3241

Singularity Theory Approach for Calculating the Runaway Boundaries of Heterogeneous Reactor Models Vaishali K. Patil, Sudhakar Subramanian, and Vemuri Balakotaiah* Department of Chemical Engineering, University of Houston, Houston, Texas 77204-4792

This work presents a systematic procedure based on singularity theory for determining the runaway boundaries of catalytic reactors. The procedure is illustrated using a well-mixed cooled heterogeneous reactor (continuous stirred tank reactor) model. Initially, a comprehensive steadystate analysis is presented and the different possible bifurcation diagrams of the solid phase temperature versus residence time are evaluated using singularity theory and the continuation technique. These results are then used to determine the different possible shapes of runaway boundaries, and a methodology is presented to calculate them. It is observed that the particle Lewis number, Lep, and the particle Damko¨hler number, Dap, are the two important parameters that determine the conditions under which the pseudohomogeneous model predictions break down. Depending on the values of these parameters, the runaway boundary is either given completely by the boundary limit set (corresponding to ignition of catalyst particles at zero residence time) or completely by the isola variety (corresponding to the coalescence of low- and high-temperature solution branches) or partially by the isola and the hysteresis (corresponding to the coalsecence of an ignition and extinction point) varieties. Finally, the influence of the reactor model type, various parameters, and time dependent perturbations on the runaway boundary is discussed. The results of this work have industrial significance because of the increasing trend of using highly active catalysts. In such cases, the interphase heat and mass transfer and intraparticle diffusion are the rate limiting steps and heterogeneous models are essential for accurate prediction of the runaway boundary. For example, unlike the pseudohomogeneous models which predict that reactor operation can be made safe by using sufficiently high cooling, heterogeneous models predict that beyond a certain critical point (boundary limit or isola variety) no amount of cooling can bring the reactor from runaway to safe region. In such cases, runaway conditions are solely determined by the transport processes occurring at the particle level. 1. Introduction Mathematical models describing the steady-state and dynamic behavior of heterogeneous reactors are obtained by writing the continuity, momentum, species, and energy balance equations and combining these with the various constitutive equations for the rate processes (for inter- and intraphase heat and mass transfer and chemical reactions). These resulting equations are usually highly nonlinear due to the strong coupling between the transport processes and the chemical kinetics and contain several parameters. It is wellknown that these nonlinear models exhibit multiple steady states, periodic states in space (and or time), and regular (or irregular) spatiotemporal patterns. Most of the earlier work on the steady-state and dynamic behavior of reactors either was focused on very simple models of these systems or was confined to numerical simulation of the governing equations. However, given the large dimension of the parameter space of these systems, a comprehensive parametric study in the multidimensional space to determine all the different possible qualitative features was unavailable. More recently, with the advances in the analysis of nonlinear systems, a more systematic approach for handling these problems has become possible. Singularity theory and the center manifold theory (Golubitsky and Schaeffer, 1985; Wiggins, 1990) are being increasingly used to * Author to whom correspondence should be addressed. E-mail address: [email protected]. S0888-5885(96)00607-0 CCC: $14.00

analyze and classify the steady-state and dynamic behavior of reacting systems. Initial applications of the singularity theory and the theory of dynamical systems was mainly focused on lumped and pseudohomogeneous models of reacting systems. The steady-state analysis of a continuous stirred tank reactor (CSTR) model with single and multiple reactions has been analyzed extensively by Balakotaiah and Luss (1983, 1984, 1986). They also provided a systematic procedure for the steady-state classification of various lumped parameter systems. Nielsen and Villadsen (1983, 1985) analyzed the steadystate multiplicity features of porous catalyst pellets for different geometries. Homogeneous reactor models described in one spatial dimension (two-point boundary-value problems) were then analyzed. The shooting method (Keller, 1972; Kubicek and Marek, 1983) combined with singularity theory has been used to study the steady-state multiplicity features. Witmer et al. (1986) presented an elegant methodology for calculating the multiplicity features for distributed problems described in one spatial dimension. This method has been applied for several reacting systems. Lovo and Balakotaiah (1992) presented examples illustrating the classification of the steady-state multiplicity behavior of convection-reaction models using the shooting method. More recently, Subramanian and Balakotaiah (1996a) presented a comprehensive study of the steady-state and dynamic classifications for several homogeneous distributed reactor models in one spatial dimension. At present, there © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3231

are very few studies that give a comprehensive analysis of the steady-state or dynamic behavior of heterogeneous (catalytic) reactor models. Once a reactor model is selected, one can apply the steady-state and dynamic singularity theory to determine the different operating states along with their stability as well as to classify the parameter space into different regions in each of which a different type of behavior exists. When this information is available, one can address practical aspects of reactor stability, such as determining the boundary between safe and runaway regions of operation. Runaway (or parametric sensitivity) of a chemical reactor describes a situation in which a small change in an operating variable such as feed temperature (or conversion) induces a large change in the temperature profile of the reactor. Runaway prevails in chemical reactors as a result of the exponential dependence of the reaction rate on the temperature and linear dependence of heat transfer or removal rates on the temperature. In other words, a small perturbation in the reaction temperature can lead to high rates of reaction due to the Arrhenius dependence of the reaction rate constant on the temperature. The increased reaction rate results in more heat generation which further increases the temperature. This feedback mechanism results in reactor runaway when the rate of heat generation far exceeds the rate of heat removal. Poor reaction yield and selectivity, rapid catalyst deactivation, and personnel safety hazards are some of the problems that accompany runaway. A considerable amount of work has been reported in the literature on reactor runaway during the past 30 years. However, most of the work has focused on homogeneous reactor models. Barkelew (1959), van Welsenaere and Froment (1970), Morbidelli and Varma (1982), and Balakotaiah (1989) have presented runaway criteria for homogeneous reactor models. Only a few studies have analyzed the behavior of heterogeneous reactor models (Morbidelli and Varma, 1986). Balakotaiah and Luss (1991) derived a simple and explicit criterion for predicting runaway in a catalytic reactor which accounts for interphase heat transfer resistance and intraparticle diffusion. Parametric sensitivity behavior of a tubular reactor with cocurrent external cooling was studied by Bauman et al. (1990). The review article by Balakotaiah and Kodra (1992) summerized the stability criteria for various types of reactors. Balakotaiah et al. (1995) applied the singularity theory to determine the runaway boundaries for different homogeneous reactor configurations with single and multiple reactions. More recently, Subramanian and Balakotaiah (1996b) presented a comprehensive steadystate and dynamic classification of a well-mixed, adiabatic, heterogeneous model using the singularity theory. The main purpose of this article is to examine the impact of the interphase heat and mass transfer resistances and cooling on the runaway boundary of a catalytic reactor. We develop a model for a well-mixed, cooled, heterogeneous reactor that accounts for the interphase heat and mass transfer resistances and intraparticle diffusion. First, we present a comprehensive steady-state classification of the model and determine the different possible types of steady-state bifurcation diagrams that exist along with the phase diagrams that classify the steady-state behavior of the system. This result is then used to determine the runaway boundaries. Our results show that the shape of the runaway boundary and its defining conditions depend strongly

on the particle Lewis number (Lep) and the particle Damko¨hler number (Dap). Unlike the homogeneous case in which the runaway boundary is determined by the balance between heat generation and heat removal rates in the reactor, the boundary for the heterogeneous case may be determined by the boundary limit set (corresponding to the ignition of catalyst particles at zero residence time) or the isola variety (corresponding to the appearance or disappearance of isolated solution branches). The practical implications of these results are that runaway may occur in a catalytic reactor at much lower values of the parameter B (dimensionless heat of reaction) than in the homogeneous case and runaway conditions may be determined by the processes occurring at the particle level rather than the reactor level (as predicted by pseudohomogeneous models). We present the model equations in the next section. We review the steady-state behavior and runaway limits of the homogeneous reactor model in section 3. In section 4, we present a brief analysis of the steady-state behavior of the heterogeneous model using singularity theory. A procedure for determining the runaway boundaries is presented in section 5. We then discuss the influence of time dependent perturbations on the runaway boundaries and the calculation of these boundaries for other heterogeneous reactor models. In the last section, we interpret the practical relevance of the results of this work. 2. Mathematical Model In this section, we develop a mathematical model for the well-mixed, cooled, heterogeneous reactor (twophase CSTR). This model is a simplified case of a fluidized bed or a slurry reactor in which the solid particles are well-mixed and suspended in the fluid. In formulating the model, we make the following assumptions: (i) a single irreversible first-order reaction with the rate constant has an Arrhenius dependence on temperature; (ii) the interphase heat and mass transfer coefficients remain constant; (iii) the relative volumes occupied by the fluid and solid phases remain constant at and 1 - , respectively; (iv) the solid phase has constant activity and is confined to the reactor, while the fluid phase enters and exits at a constant flow rate; and (v) the particles are small enough so that no temperature gradients exist within the solid phase. We consider both the solid and fluid phases separately and write the species and energy balance equations for both. The governing equations in dimensionless form can be obtained as

(1 - )

( )

γθs xs - xf dxs )+ (1 - xs)η exp dt Dapm γ + θs

£ σ(1 - )

xf xs - xf dxf )+ dt Da Dapm

dθf -θf θs - θf ) + - R(θf - θc) dt Da Daph

( )

dθs γθs θs - θf )+ B(1 - xs)η exp dt Daph γ + θs

(1)

Some of the important dimensionless numbers and parameters that appear in the above model are as follows:

3232 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

xf )

Cfo - Cf Cfo θs ) γ

Ts - Tfo Tfo R)

B)

Cfo - Cs Cfo

xs )

Tf - Tfo Tfo

3. Review of the Steady-State Classification and Runaway Boundaries for the Pseudohomogeneous Model

θf ) γ

Tc - Tfo Tfo

The limiting case of negligible particle heat and mass transfer resistances (Dapm f 0, Daph f 0) reduces the heterogeneous model to a pseudohomogeneous model. The governing equations are obtained from eq 1 as

θc ) γ

UA FfCpf(1 - )k(Tfo)

dx x γθ )+ (1 - x) exp η dt Da γ+θ

(

tc (1 - )k(Tfo)VR ) tr q E γ) RTfo

γ(-∆H)Cfo FfCpfTfo

Da )

Le*

tm (1 - )k(Tfo) ) tr kcav th (1 - )k(Tfo)FfCpf Daph ) ) tr hav

Dapm )

Lep )

Dapm h ) Daph FfCpfkc

σ)

FsCps FfCpf

£ ) Le + σ(1 - )(Le - 1) (Mcp)w Le ) 1 + VR[Ffcpf + (1 - )Fscps] η)

3φp coth(3φp) - 1 2

3φp

φpo2 )

k(Tfo) Vp De Sp

t)

t′ tr

k(Ts) ) k(Tfo) exp

tc ) tm )

( ) ( )

φp2 ) φpo2 exp

VR q

1 kcav

γθs γ + θs

(

x + (1 - x) exp(θ) ) 0 Da

(4a)

θ + B(1 - x) exp(θ) - Rθ ) 0 Da

(4b)

-

-

The reactor Damko¨hler number (Da) is the ratio of residence time (tc) to that of the characteristic reaction time (tr). The parameter γ is the dimensionless activation energy of the reaction and is a measure of the sensitivity of the reaction rate to changes in temperature. The parameter B is the dimensionless adiabatic temperature rise and represents the maximum temperature attainable by the fluid under adiabatic conditions. The parameters Dapm and Daph represent the particle mass and heat Damko¨hler numbers, respectively. The parameter Dapm is the ratio of the time constant for external mass transfer (tm) to that of the reaction, while the parameter Daph is the ratio of the time constant for solid-fluid heat transfer (th) to that of the reaction. A small value for Dapm (Daph) implies negligible mass (heat) transfer resistance and vice versa. The particle Lewis number, Lep, is the ratio of Dapm to Daph, or, equivalently, the time constant for external mass to heat transfer. The parameter σ denotes the ratio of the volumetric heat capacities of the solid to that of the fluid, φp is the normalized Thiele modulus, and the reactor Lewis number, Le, is the ratio of the heat capacity of the reactor to that of the reactor contents. The parameter η is the internal catalyst effectiveness factor and is unity if no concentration gradients exist within the particles. Finally, we note that while the general (transient) model contains eleven parameters, the steady-state model contains only eight parameters.

)

Here, Le* ) £ + σ(1 - ). Note that since there is no distinction between the fluid and solid phases, there is a single temperature (θ) and conversion (x). This model reduces to the homogeneous CSTR model in the limit of η ) 1 and ) 1. In this section, we present a brief review of the homogeneous model to illustrate the use of singularity theory to determine the phase diagrams and runaway boundaries. The above model has been extensively analyzed in the literature, and more details can be found elsewhere (e.g., Balakotaiah and Luss (1984)). Since, our aim in this section is to present the main ideas involved, we further simplify the model (given by eq 3) by assuming (i) infinite activation energy (γ f ∞) and (ii) the coolant temperature to be the same as the feed temperature (θc ) 0). The equations at steady state can now be written as

tr )

(2)

(3a)

θ γθ dθ )+ B(1 - x) exp η - R(θ - θc) dt Da γ+θ (3b)

γθs γ + θs

1 (1 - )k(Tfo) Ffcpf th ) hav

)

Equations 4a and 4b can be combined and rewritten as

F(θ,Da,B,R) ) 1 + R + [B - θ(1 + RDa)] exp(θ) ) 0 (5) -θ Da

(

)

The above equation can be shown to have a maximum of three solutions. The details of obtaining the steadystate classification for equations such as in eq 5 (dependent on a single state variable) can be found in Balakotaiah and Luss (1984). The steady-state classification classifies the parameter space into different regions such that the steady-state bifurcation diagram (plot of state variables such as temperature versus bifurcation parameters such as Da) in each region is different. For the above model, the steady-state classification diagram is obtained by computing the following varieties. (i) Hysteresis Variety (H). This variety is obtained by solving the following three equations:

F)

∂F ∂2F ) )0 ∂θ ∂θ2

(6)

This variety divides the parameter space into two regions, one in which there is a unique steady state for all Da values and one where the bifurcation diagram contains a hysteresis loop. Crossing this variety, an ignition and an extinction point either appear or disappear from the steady-state solution branch. From eqs 5 and 6, the hysteresis variety can be obtained as


e2 (7) x ) 0.5, θ ) 2, Da ) exp(-2), R ) (B - 4) 4 (ii) Isola Variety (I). This variety is obtained by solving the following three equations:

F)

∂F ∂F ) )0 ∂θ ∂Da

(8)

Crossing this variety in the parameter space, an isolated solution branch either appears (or disappears) or two solution branches either join or separate in the steadystate bifurcation diagram. From eqs 5 and 8, the isola variety can be written in a parameteric form as 2

B)

(1 - x) 1 1 exp , R) 2 2 1 x x (1 - x) x

(

)

0 2.0 and decreases from (8,0) to (4Lep,e-2) for Lep < 2.0. For Lep ) 2.0, the pitchfork is independent of Dap and is just a horizontal line. The locus in the (R,Dap) is however not monotonic for all Lep values. For 0 < Lep < 0.5, the pitchfork locus intersects the R ) 0 axis at a finite Dap value (Dapc). This implies that, for

∂F ∂F )0 ) ∂θs ∂Da

at Da ) 0

(22)

Solving eq 22 along with eq 13, we can show that this locus can be expressed as

θs ) B ) 1/(1 - Lep) Dap )

(

) (

)

Lep -1 exp 1 - Lep 1 - Lep

(23)

As expected, this locus is independent of the cooling parameter, R. This locus gives us the critical parameter values (B and Dap) above (or below) which there is either no isola present or one of the branches of the isola goes into the nonfeasible range (this will get clearer when the classification plots are shown). Hence, for a fixed


value of Lep, we can calculate the corresponding critical B and Dap values from eq 23. It is observed that this locus exists only for Lep < 1.0. Hence, for Lep > 1.0, this higher codimension singularity does not influence the qualitative nature of the classification diagram. The determination of the critical parameter values obtained from these higher order (codimension 2) singularities will be explained when we discuss the steady-state classification plots. We now present the different steady-state phase plots in the (R,B) plane for fixed values of Lep and Dap. It should be pointed out that the phase diagrams can be plotted in the (B, Dap) or any other plane. From the higher codimension singularities discussed in this section, we can identify the following three different Lep ranges: (i) 0 < Lep < 0.5, (ii) 0.5 < Lep < 1.0, and (iii) Lep > 1.0. Within each region we now determine the Dap values in each of which the steady-state classification diagrams are qualitatively different. Case 1: 0 < Lep < 0.5. We will first identify the different ranges of Dap values in each of which the steady-state classification plot is qualitatively different. To obtain this information, the higher order singularities along with the boundary limit sets are plotted in the (B, Dap) plane (not shown here). Through examination of this plot, the following different Dap ranges are obtained: (i) Dap3 < Dap < exp(-2), where no pitchfork exists, no isola variety exists, and the hysteresis variety is to the right of the boundary limit sets (Figure 3a); (ii) Dap2 < Da < Dap3, where no pitchfork exists, isola variety exists, and the hysteresis variety lies to the right of the boundary limit sets (Figure 3b); (iii) Dap1 < Da < Dap2, where no pitchfork exists, isola variety exists, and the hysteresis variety lies partially between both the boundary limit sets (Figure 3c); (iv) Dap0 < Da < Dap1, where both pitchfork and the isola variety exist and the hysteresis variety lies partially between both boundary limit sets (Figure 3d); and (v) 0 < Da < Dap0, which is the same as in case iv but the hysteresis variety now intersects the R ) 0 axis to the left of both boundary limit sets (Figure 3e). The critical Dap values mentioned above are obtained as follows. The Dap3 value is obtained from eq 23. The hysteresis variety evaluated at R ) 0 intersects the boundary limit set at two points, Dap2 and Dap0. The critical value Dap1 corresponds to the Dap value where the pitchfork locus intersects the R ) 0 axis. The five different classification diagrams corresponding to the above Dap ranges are summarized in Figure 3. The transitions in the bifurcation diagrams as we move from one region to the other are quite easy to comprehend and are not discussed here. Case 2: 0.5 < Lep < 1. As in the earlier case, we first determine the critical Dap values using the higher codimension singularities. The pitchfork locus, unlike the previous case, exists for the entire range of Dap values (0 < Dap < e-2). The Dap ranges for this case can be summarized as follows: (i) Dap5 < Dap < exp(-2), where the isola variety has two branches and the pitchfork point exists and is to the left of the two boundary limit sets (Figure 4a); (ii) Dap4 < Da < Dap5, where the isola variety has two branches and the pitchfork point exists and is between the two boundary limit sets (Figure 4b); and (iii) 0 < Da < Dap4, where the isola variety has a single branch and the pitchfork point is between the two boundary limit sets (Figure 4c). The critical Dap values mentioned above are obtained as follows. In the (B, Dap) plane, for Dap > Dap5, the B value at the pitchfork point lies to the left of the two

Figure 4. Schematic classification diagrams of the steady-state behavior in the (R,B) plane for 0.5 < Lep < 1.0 and (a) Dap5 < Dap < exp(-2), (b) Dap4 < Dap < Dap5, and (c) 0 < Da < Dap4.

boundary limit sets, while, for Dap < Dap5, the pitchfork is between the two boundary limit sets. Thus, at Dap ) Dap5, the pitchfork coordinates fall on the extinction boundary limit set. The value of Dap4 is obtained from eq 23. The phase diagrams for these three cases are summarized in Figure 4. Case 3: Lep > 1. For this range of Lep values there is a single classification diagram for the entire Dap range. When we consider the higher codimension singularities, we observe that there are no critical Dap values. The classification diagram is qualitatively similar to that in Figure 4b. The numerically computed bifurcation diagrams in each of the regions of Figures 3 and 4 are shown in Figure 5 for one set of parameter values in each region. The above analysis on the steady-state behavior will now be used to determine the runaway boundaries. 5. Determination of Runaway Boundaries For a heterogeneous reactor model, excursions to the high-temperature branch from the low-temperature branch due to small perturbations in the operating variables can occur in either of the following two ways: (i) The reaction path connecting xf ) 0, θs ) 0 and xf ) 1, θs ) θc contains ignition points, or (ii) the heat and mass transfer resistances between the solid and fluid phases are such that the particle temperature (θs) exists


reaction path connecting xf ) 0, θs ) 0 to xf ) 1, θs ) θc has no ignition points and in addition, along this reaction path, the particle equation has a low-temperature steady-state solution. As shown elsewhere (Balakotaiah and Kodra, 1992; Balakotaiah et al., 1995) this definition of the runaway boundary is consistent with the pseudohomogeneous models and coincides with the parametric sensitivity locus defined in the literature. Moreover, as will be shown later, this locus is also independent of the reactor model; i.e., it is almost identical for all heterogeneous reactor models. Most of the analysis in the literature to determine the runaway boundaries have used the pseudohomogeneous models. For the heterogeneous reactor models, as defined above, the existence of ignition points arising due to the particle equations also needs to be taken into account. We first consider the different ranges of Lep and Dap values, as discussed in section 4, to construct the runaway boundaries in each region. Finally, we integrate the results to present a systematic procedure for determining the runaway boundaries. Case 1: 0 < Lep < 0.5. For this range of Lep values, depending on the Dap values, five different phase diagrams were presented in section 4. However, for the determination of runaway boundaries, we can reduce that number on the basis of inspection of the classification plots (Figure 3), to analysis of three ranges of Dap values: (a) Dap3 < Dap < e-2. A schematic of the classification diagram for this case was shown in Figure 4. From this plot, we observe that, in regions i and ii, the lowtemperature branch is stable for small perturbations. However, in region iii an ignition point exists and only a high-temperature state exists for low Da values. Similarly, in region iv, the temperature for low Da values is high (refer to Figure 5). Hence, the runaway boundary for this case is given by the ignition boundary limit set which separates regions i and ii (safe regions) from regions iii and iv (unsafe regions). This boundary can be expressed in a parametric form as

B)

Lep xs(1 - xs)

0 < xs < 0.5

(24)

where

Dap ) Dap3 )

Figure 5. Numerically computed bifurcation diagrams corresponding to regions i-xi in Figures 3 and 4. The other parameter values are Lep ) 1.0 and Dap ) 0.085 for cases i-vi and viii-xi. For case vii, Lep ) 0.2 and Dap ) 0.01.

in an ignited state even though the temperature of the fluid (θf) is low. Hence, we define the operation of a heterogeneous reactor to be insensitive to small perturbations if the

(

(

xs -1 exp 1 - xs 1 - xs

) (

)

(25)

)

Lep -1 exp 1 - Lep 1 - Lep

(26)

It should be noted that this boundary is independent of the cooling parameter, R, and hence the runaway in the (R/B, B) plane is a straight line parallel to the R axis. Thus, for this range of Dap values, obtaining the runaway locus involves calculating the critical B values using eqs 24 and 25. The runaway locus is shown in Figure 6a for two different Dap values and Lep ) 0.2. The physical interpretation for this case is as follows: for B < Bcritical, the reactor is in the safe region irrespective of the degree of cooling. For B > Bcritical, the reactor is in the region of runaway and no amount of cooling can bring the reactor from a runaway state to a safe state. This is due to the fact that, for this range of Dap and Lep values, the runaway is determined by the transport processes occurring at the particle level. (b) Dap1 < Dap < Dap3. The classification diagram corresponding to this case is shown in Figure 3. Here,


Figure 7. Numerically computed runaway locus for the heterogeneous model for Lep ) 0.2 and Dap values in the range 0 < Dap < Dap1.

Figure 6. Numerically computed runaway locus for the heterogeneous model for Lep ) 0.2 and (a) Dap3 < Dap < exp(-2) and (b) Dap1 < Dap < Dap3.

due to the presence of the isola and hysteresis varieties, region v, which has ignition points, exists. Hence, the runaway locus which separates regions i and ii from the unsafe regions is given completely by the isola variety (Figure 3). The hysteresis and boundary limit varieties do not influence the runaway locus. The runaway boundary can be represented in a parametric form by eq 17. The value of Dap1 corresponds to the intersection of the pitchfork locus with the R ) 0 axis. This critical value can be shown to be given by the expression

Dap1 )

(

)

Lep exp(-2) 1 - Lep

(27)

The runaway locus for this range of Dap values is shown in Figure 6b for three different Dap values and Lep ) 0.2. It can be seen that the runaway boundary has moved toward lower B values. (c) 0 < Dap < Dap1. The classification plots corresponding to this range of Dap values are shown in Figure 3. Due to the presence of regions vii and viii, the corresponding bifurcation diagrams (Figure 5) contain ignition points and the runaway locus is now given by the hysteresis variety for lower B values and by the isola variety for B values beyond the pitchfork point. The runaway locus is given as

B)

4Lep(1 + R(e-2 - Dap))

, Lep + e2Dap(1 - Lep + R(e-2 - Dap)) for B < B* (28)

and by eq 17 with 1 < θs< 2 for B > B*. The value of B* is obtained from eqs 20c and 20e.

Figure 8. Comparison of the numerically computed runaway boundaries of the heterogeneous model for Lep ) 1.0 and different Dap values with the pseudohomogeneous model.

The runaway locus for this range of Dap values is shown in Figure 7 for varying Dap values and Lep ) 0.2. Case 2: Lep > 0.5. The schematic classification plots for this range of Lep values is shown in Figure 4. Hence, the only new feature is that the isola variety, depending on the Dap value, can have two branches. However, by our definition of runaway locus, we observe that the new region (ix) is safe for all practical purposes. The runaway boundary is once again given partially by the hysteresis and isola varieties. The expressions are the same as those given in eqs 17 and 28. Figures 8 and 9 present a few runaway loci for different values of Dap and Lep values of 1.0 and 2.0, respectively. It is interesting to note that, for Lep > 1, the pseudohomogeneous model predicts runaway in the safe regions (for low R values). On the basis of the above discussion, the runaway boundary can be evaluated as follows: (i) determine the values of Lep and Dap using eq 2; (ii) if Lep > 0.5, the runaway boundary is calculated using eqs 17 and 28; (iii) if Lep < 0.5, the values of Dap3 and Dap1 are determined from eqs 26 and 27. The runaway boundary can then be determined using either eqs 24 and 25 or eq 17 or eqs 17 and 28. We now discuss how the runaway boundaries for the pseudohomogeneous and heterogeneous models compare with each other. More importantly, it needs to be determined under what conditions the predictions of the pseudohomogeneous model break down. Figures 8 and 9 show the runaway boundaries for the pseudohomogeneous model and the heterogeneous model. One


Figure 9. Comparison of the numerically computed runaway locus of the heterogeneous model for Lep ) 2.0 and different Dap values with the pseudohomogeneous model. Table 1. Comparison of the Critical Lep Values between the Simplified and the General Models γ and φpo values

range of Lep values

γ f ∞, φpo ) 0.0

0 < Lep < 0.5 0.5 < Lep < 1.0 Lep > 1.0 0 < Lep < 0.435 0.435 < Lep < 0.682 Lep > 0.682 0 < Lep < 0.385 0.385 < Lep < 0.675 Lep > 0.675 0 < Lep < 0.305 0.305 < Lep < 0.47 Lep > 0.47

γ ) 25, φpo ) 0.1 γ ) 25, φpo ) 1.0 γ ) 15, φpo ) 1.0

critical Dap value exp(-2.0) 0.108 046 0.110 307 0.078 353

obvious difference between the two reactor models can be readily observed. The pseudohomogeneous model predicts that, for significantly high cooling, the reactor always exists in a safe region. On the other hand, the heterogeneous model predicts that, for B > Bcritical (a critical value), no amount of cooling can bring the reactor operation to a safe region. As the Dap value increases, the deviation between the predictions of the two models also increases. The vertical asymptote (B ) Bcritical) arises because, under certain conditions, the solid particles can ignite (at Da ) 0 or zero residence time), while the fluid temperature can still be low. Thus, one should exercise caution when using the pseudohomogeneous model to predict runaway boundaries for catalytic reactors. 6. Discussion The analysis presented so far assumed infinite activation energy, neglected the presence of intraparticle resistances, and did not take into account the existence of time dependent periodic solutions. In this section, we analyze how the parameters as well as the presence of time dependent solutions and model type influence the runaway boundaries. Influence of Activation Energy and Intraparticle Resistances on the Runaway Boundaries. In order to determine the influence of these parameters on the runaway boundaries, a detailed steady-state classification of the reactor model was carried out. The main conclusion is that the qualitative behavior of the steady-state classification remains unaffected by varying these parameters. However, we also observed that the values of the critical parameters (Lep, Dap) were different depending on the values of the activation energy, γ, and the Thiele modulus, φpo. Table 1 gives

Figure 10. Numerically computed runaway boundaries illustrating the influence of intraparticle diffusional resistance and the activation energy.

the critical Lep ranges along with the critical Dap value for different values of γ and φpo (beyond the critical Dap value, there always exists a unique high-temperature solution). An assumption of unity for the effectiveness factor implies that the concentration of the reactant inside the solid phase is uniform (φpo ) 0). Hence, the heat released for finite φpo values will always be less than that for the simplified case of uniform concentration. Thus, the runaway boundary in the (R/B, B) plane for the real case will shift toward the right, implying that the simplified model predicts a part of the safe region to be unsafe. A similar reasoning shows that the assumption of infinite activation energy results in a conservative prediction of the runaway boundary. We note that in the strong pore diffusional limitation case (φ . 1 so that ηφ ) 1) the analysis presented in sections 4 and 5 can be repeated with

θs* )

θs , Da* ) Da/φpo, Dap* ) Dap/φpo, 2 B* ) B/2, R* ) φpoR (29)

Thus, the runaway boundary for the two limiting cases of no pore diffusional limitation and strong intraparticle resistance can be determined analytically. Numerical calculations are needed only for the intermediate cases of finite γ and φpo. We show in Figure 10 some numerically computed runaway boundaries for finite values of γ and φpo > 0, along with the simplified cases. As mentioned earlier, we observe that, when the value of γ is decreased or the value of φpo is increased, the runaway boundaries shift to the right (higher B values). Influence of Transient Solutions on the Runaway Boundaries. Time dependent solutions arise in the system when the steady-state solution branch loses stability through a Hopf bifurcation (a pair of purely imaginary eigenvalues cross the imaginary axis at a nonzero speed). In many cases, it is undesirable to have oscillatory solutions arising in the system. A detailed dynamic classification of the reactor model must be


performed to determine the conditions leading to oscillatory solutions as well as to determine the influence on the runaway boundaries. A detailed analysis of this aspect, for several homogeneous reactor models and a well-mixed heterogeneous reactor model, can be found in Subramanian and Balakotaiah (1996a,b). The main results of the study were as follows: (i) In the feasible parameter space, high-temperature oscillations are more likely to occur. Since these occur in the unsafe regions, these do not influence the runaway boundaries obtained in section 5. (ii) In situations where lowtemperature oscillations occur, the residence times are very high. Moreover, the amplitudes of these oscillations are not high enough to cause runaway in the reactor. Hence, it is concluded that the runaway boundaries predicted using the steady-state analysis represent the stability boundary for stationary as well as time dependent perturbations. Influence of the Model Type on the Runaway Boundaries. The runaway boundaries presented here were obtained using the cooled, well-mixed model. It is important to determine how varying the reactor configuration alters the runaway boundaries. From the results presented in this section, we observe that the runaway boundaries for heterogeneous models consist of two asymptotes: (i) the adiabatic and (ii) the particle ignition asymptotes. The particle ignition asymptote is given by the value of B corresponding to the ignition boundary limit set. It is also seen from the definition of the boundary limit sets that it corresponds to particle ignition at zero residence time (or inlet fluid conditions). In other words, it is independent of the specific type of reactor. The adiabatic asymptote for the pseudohomogeneous and heterogeneous models (irrespective of the reactor configuration) are very close to each other. This result, for different pseudohomogeneous reactor models, has been shown by Balakotaiah et al. (1995). Hence, irrespective of the model type, a reasonably accurate runaway boundary can be obtained using the wellmixed, heterogeneous reactor model. In most practical situations, the uncertainity in the kinetic and transport parameters is much larger than the differences between various reactor models. Thus, the runaway boundaries presented here can be considered to be “universal” and applied to all catalytic reactors such as multitube fixed bed reactors, fluidized bed reactors, slurry reactors, and heterogeneous batch reactors. 7. Conclusions The main objective of this work was to determine the runaway boundaries for heterogeneous reactor models and compare the results with that of the pseudohomogeneous reactor models. A well-mixed, cooled, heterogeneous reactor model is developed and analyzed. First, a detailed steady-state classification of the reactor model is performed. The results indicate that depending on the parameter values we obtain eight different classification diagrams and twelve bifurcation diagrams. These results are then used to outline a procedure for calculating the runaway boundaries. It is observed that the particle Lewis number Lep (ratio of time constant for mass to heat transfer) and the particle Damko¨hler number Dap values determine how the runaway boundaries are evaluated. Depending upon the values of these parameters, the runaway boundary is either given by (i) entirely the boundary limit set, (ii) entirely by the isola variety, or (iii) partially by the isola and the hysteresis

varieties. For Lep > 0.5, irrespective of the Dap value, the runaway locus is given by the hysteresis and the isola varieties. For Lep < 0.5, the value of Dap determines how the runaway boundary is defined. When the Dap value is small, it is defined by the hysteresis and the isola variety. As the value of Dap is increased, the runaway boundary is given entirely by the isola variety, and, for large values of Dap (close to e-2), the runaway boundary is given by the ignition branch of the boundary limit set. This is different from the pseudohomogeneous case where the runaway boundary is entirely given by the isola and the hysteresis varieties. The influence of finite activation energy and intraparticle resistances is then discussed. It is reasoned that the runaway boundaries predicted by the simplified model are conservative and, for practical purposes, the runaway boundaries calculated using the simplified model are accurate enough. The time dependent perturbation are also shown to have negligible influence in the determination of the runaway boundaries. Even though the runaway boundaries are obtained using the simplified model, it is concluded that changing the reactor configuration will not alter the runaway boundaries. The runaway boundaries presented in this paper can be used to determine the safe and unsafe regions of operation for most reactor operations. For situations which warrant exact prediction of runaway boundaries, these simplified results can be used as a good starting point for performing detailed computations using specific kinetics and reactor models. Acknowledgment This work was partially supported by a grant from the Robert A. Welch Foundation. Nomenclature A: heat transfer area (m2) av: interfacial area per unit reactor volume (m-1) B: dimensionless adiabatic temperature rise CA: concentration (mol‚m-3) cp: specific heat (J‚kg-1‚K-1) Daph: particle Damko¨hler number for heat Dapm: particle Damko¨hler number for mass dp: diameter of particle (m) E: activation energy for the reaction (J‚mol-1) -∆H: heat of reaction per mole of extent (J‚mol-1) h: fluid-particle heat transfer coefficient (W‚m-2‚K-1) k(T): reaction rate constant (s-1) kc: mass transfer coefficient (m‚s-1) Le: Lewis number p*: vector of parameters q: volumetric flow rate of fluid (m3‚s-1) r: rate of reaction per unit volume of the solid phase (mol‚m-3‚s-1) Sp: external surface area of catalyst particle (m2) T: temperature (K) t: dimensionless time U: overall heat transfer coefficient (W‚m-2‚K-1) Vp: volume of catalyst particle (m3) VR: volume of reactor (m3) Greek Letters : porosity/fractional volume occupied by fluid φ2: Thiele modulus γ: dimensionless activation energy λ: thermal diffusivity (m2‚s-1) µ: viscosity of fluid (kg‚m-1‚s-1) θ: dimensionless temperature Ff: density of the fluid (kg‚m-3)

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3241 σ: ratio of heat capacities of solid to fluid Subscripts and Superscripts c: cooling critical: critical condition eff: effective f: fluid phase h: heat m: mass o: inlet conditions p: particle R: reactor s: solid phase w: reactor walls

Literature Cited Balakotaiah, V., Simple Runaway Criteria for Cooled Reactors. AIChE J. 1989, 35, 1039. Balakotaiah, V.; Luss, D. Multiplicity Features of Reacting Systems. Chem. Eng. Sci. 1983, 38, 1709. Balakotaiah, V.; Luss, D. Global Analysis of the Multiplicity Features of Multiple Reaction Lumped Parameter Systems. Chem. Eng. Sci. 1984, 39, 865. Balakotaiah, V.; Luss, D. Steady-State Multiplicity Features of Lumped Parameter Chemically Reacting Systems. In Dynamics of Nonlinear Systems; Hlavacek, V., Ed.; Gordon and Breach: New York, 1986. Balakotaiah, V.; Luss, D. Explicit Runaway Criterion for Catalytic Reactors with Transport Limitations. AIChE J. 1991, 37, 1780. Balakotaiah, V.; Kodra, D. Stability Criteria for Chemical Reactors. In Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes (Dycord+ 92 preprints), College Park, MD, Balchen, J. G., Ed.; Pergamon: Oxford, U.K., 1992. Balakotaiah, V.; Kodra, D.; Nguyen, D. Runaway Limits for Homogeneous and Catalytic Reactors. Chem. Eng. Sci. 1995, 50, 1149. Barkelew, C. H. Stability of Chemical Reactors. Chem. Eng. Prog., Symp. Ser. 1959, 55, 38. Bauman, E.; Varma, A.; Lorusso, J.; Dente, M.; Morbidelli, M. Parametric Sensitivity in Tubular Reactors with Co-current External Cooling. Chem. Eng. Sci. 1990, 45, 1301.

Golubitsky, M.; Schaeffer, D. G. Singularities and Groups in Bifurcation Theory; Springer-Verlag: New York, 1985; Vol. I. Keller, H. B. Numerical Methods for Two Point Boundary Value Problems; Blaisdell: New York, 1972. Kubicek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer: New York, 1983. Lovo, M.; Balakotaiah, V. Multiplicity Features of Adiabatic Autothermal Reactors. AIChE J. 1992, 38, 101. Morbidelli, M.; Varma, A. Parametric Sensitivity and Runaway in Tubular Reactors. AIChE J. 1982, 28, 705. Morbidelli, M.; Varma, A. Parametric Sensitivity in Fixed-Bed Reactors: the Role of Interparticle Transfer Resistances. AIChE J. 1986, 32, 297. Nielsen, P. H.; Villadsen, J. Absorption with Exothermic Reaction in a Falling Film Column. Chem. Eng. Sci. 1983, 38, 1439. Nielsen, P. H.; Villadsen, J. An Analysis of the Multiplicity Pattern of Models for Simultaneous Diffusion, Chemical Reaction and Absorption. Chem. Eng. Sci. 1985, 38, 1439. Subramanian, S.; Balakotaiah, V. Classification of Steady-State and Dynamic Behavior of Distributed Reactor Models. Chem. Eng. Sci. 1996a, 51, 401. Subramanian, S.; Balakotaiah, V. Classification of Steady-State and Dynamic Behavior of Well-Mixed Heterogeneous Reactor Model. Chem. Eng. Sci. 1996b, 52, 961. van Walsenaere, R. J.; Froment, G. F. Parametric Sensitivity and Runaway in Fixed-Bed Catalytic Reactors. Chem. Eng. Sci. 1970, 25, 1503. Wiggins, S. Introduction to Applied Non-linear Dynamical Systems and Chaos; Springer-Verlag: New York, 1990. Witmer, G.; Balakotaiah, V.; Luss, D. Finding Singular Points of Two-Point Boundary Value Problems. J. Comput. Phys. 1986, 65, 244.

Received for review October 2, 1996 Revised manuscript received March 31, 1997 Accepted April 12, 1997X IE960607H

X Abstract published in Advance ACS Abstracts, July 1, 1997.

Singularity Theory Approach for Calculating the Runaway Boundaries

Recommend Documents