Ind. Eng. Chem. Res. 1987, 26, 1351-1356 Experimental solvent mixture dielectric constant data of Akerlof (1932) are used to represent dielectric constant as a function of concentration D = AI AzXf3 A3X’32 A4XfS3 A6Xi4 A6X’35 03-8)
+
+
+
+
+
The temperature dependency of the dielectric constant at a constant composition is represented by In D = 2.303[AD1 + A&(T - 293.15)] (B-9) The constants AD1 and AD2 are taken from Akerlof (1932).
Literature Cited Akerlof, G. J. Am. Chem. SOC.1930,52,2353. Akerlof, G. J. Am. Chem. SOC.1932,54,4125. Beckerman, E.; Tassios, D. ACS Adv. Chem. Ser. 1976,155,3. Bixon, E.; Gurrey, R.; Tassios, D.J. Chem. Eng. Data 1979,I, 9. Bromley, L. A. J. Chem. Thermodyn. 1972,4,669. Bromley, L. A. AZChE J. 1973,19,13. Chen, C. C.; Britt, H. E.; Boston, J. F.; Evans, L. B. ACS Symp. Ser. 1980,133,61. Ciparis, J. N. Data of Salt Effects in Vapor Liquid Equilibrium; Edition of Lithuaniun Agricultural Academy: Kaunas, USSR, 1966. Covington, A. K.; Dickinson, T. Physical Chemistry of Organic Solvents Systems; Plenum: London and New York, 1973. Cruz, L. L.; Renon, H. AZChE J. 1978,24,817. Debye, P.; HOckel, E. Phys. 2.1923,24,334. Gronwall, T. H.; Lamer, V. K.; Greiff, L. J. J. Phys. Chem. 1931,35, 2245. Granwall, T. H.; Lamer, V. K.; Sandved, K. Phys. Z. 1928,29,358. Guggenheim, E.A.; Turgeon, J. C. Phil. Mag. 1935,7,585. Hala, E. Ind. Chem. Eng. Symp. Ser. 1969,32,5. Harned, H. S.; Owen, B. B. The Physical Chemistry of Electrolytic Solutions, 3rd ed.; ACS Monograph Series No. 137;Reinhold: New York, 1958;Appendix A, p 719.
1351
Marina, J. M.; Tassios, D. Znd. Eng. Chem. Process Des. Deu. 1973, 12, 67. McGlashan, M. L.; Williamson, A. G. J. Chem. Eng. Data 1976,21,
_196. _ _.
Meissner, H. P.; Kusik, C. L. AZChE J. 1972,18,294. Meissner, H. P.; Kusik, C . L.; Tester, J. W. AZChE J . 1972,18,661. Mock, B.; Evans, L. B.; Chen, C. C. AIChE J. 1986,32,1655. Perry, R. H.;Chilton, C. H. Chemical Engineers’ Handbook, 5th ed.; McGraw Hill: New York, 1973. Pitzer, K. S. J. Phys. Chem. 1973,77,268. Pitzer, K. S., Acc. of Chem. Res., 1977,10,371 Pitzer, K.S.;Mayorga, G. J. Phys. Chem. 1973, 77,2300. Pitzer, K. S.;Mayorga, G. J. S o h . Chem. 1974,3,539. Pitzer, K.S.; Kim, J. J. J . Am. Chem. SOC.1974,96,5701. Prausnitz, J. M.; Eckert, C. A.; Orye, R. V.; O’Connel, J. P. Computer Calculations for Multicomponent Vapor-Liquid Equilibria, Monograph: Prentice-Hall: Englewood Cliffs, NJ, 1967. Prausnitz, J. M.; Edwards, T. J.;Newman, J. AZChE J. 1975,21,248. Prausnitz, J. M.; Edwards, T. J.; Maurer, G.; Newman, J. AZChE J. 1978,24,966. Rastogi, A. Ph.D. Dissertation, New Jersey Institute of Technology, Newark, 1981. Rebolleda, F.; Ocon, J. An. R. SOC.Esp. Fis. Quim., Ser. E 1958,544 525. Renon, H.; Prausnitz, J. M. AZChE J. 1968,14,135. Robinson, R. A,; Stokes, R. H. Electrolyte Solutions; Butterworths: London, 1955. Rousseau, R. A.; Ashcraft, D. L.; Schoenborn, E. M. AZChE J. 1972, 17, 828. Rousseau, R. A.; Boone, J. E.; Schoenborn, E. M. ACS Adu. Chem. Ser. 1976,155,36. Rousseau, R. A.; Boone, J. E. AZChE J. 1978,24,718. Sander, B.; Fredenslund, A.; Rasmussen, P., Paper Presented at the AIChE Meeting, San Francisco, CA, Nov 1984. Skabichevskii, P. A. Russ. J. Phys. Chem. 1969,43,1432. Taniguichi, H.; Janz, G. J. J. Phys. Chem. 1957,61, 688.
Received for review February 3, 1986 Accepted March 12, 1987
Kinetics of Pyrolysis of Moroccan Oil Shale by Thermogravimetry Deepak S. Thakur* HarshawlFiltrol Partnership, Research & Development, Beachwood, Ohio 44122
H. Eric Nuttall, Jr. Department of Chemical and Nuclear Engineering, University of New Mexico, Albuquerque, New Mexico 87131
T h e kinetics of the thermal decomposition of Moroccan oil shale have been studied by isothermal and nonisothennal thermogravimetry (TG). The nonisothermal weight loss data have been analyzed by using three models, (1)Chen and Nuttall model, (2) Coats and Redfern model, and (3) Anthony and Howard model, while the isothermal TG data have been correlated by using the integral method of kinetic analysis. The combined use of nonisothermal and isothermal T G measurements has shown that the thermal decomposition of Moroccan oil shale involves two consecutive reactions with bitumen as a n intermediate. The pyrolysis rates of Moroccan shale are compared with those of Colorado and Jordan shale.
1. Introduction Morocco has moderate deposita of oil shale which can be exploited using American processing technology. To date, U.S. processes are based on the thermal treatment of the shale to decompose the insoluble organic matter called kerogen into the desired oil product. The key step in the thermal treatment is the pyrolysis stage and involves the kinetic rate of kerogen decomposition and the corresponding oil formation. Because oil shales have differing origins and are found in differing geological environments, 0888-5885/87/2626-1351$01.50/0
it is not surprising that they behave differently when subjected to pyrolysis conditions. This study was initiated to investigate the pyrolysis kinetics of Moroccan oil shale. Numerous attemps (Hubbard and Robinson, 1950; Allred 1966; Dericco and Barrick, 1956) dealing with the mechanistic and kinetic points of view have been made to understand the processes occurring during oil shale pyrolysis. Hubbard and Robinson (1950) employed an isothermal technique to obtain the pyrolysis rate data on powdered samples. The main drawback of this method is the time
0 1987 American Chemical Society
1352 Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987 Table I. Chemical Analysis of Moroccan Oil Fischer oil assay, L/ton density, g/mL chem composition of raw shale compd, wt 7' 0 .4.kerogen
c
H
S N 0 B. dolomite C. calcite D. quartz E. argiles (illite-kaolinite) F. pyrites G. FeC03 + Fe203 H. TiOz + phosphate I. other constituents
Shale 71 2.25
9.5 69.7 6.7 9.1 2.3 12.2 (by diff) 14.4 37.7 17.2 12.2 1.6 1.6 2.2 3.6
required for the sample to reach from room temperature to the reaction temperature. Braun and Rothman (1975) recognized this fact and reanalyzed the results of Hubbard and Robinson (1950) by making corrections to account for the thermal induction period. Recently, Smith and coworkers (Galan and Smith, 1983; Pan, et al., 1985) determined the intrinsic kinetics of Colorado oil pyrolysis by incorporating induction period. Nonisothermal thermogravimetric analysis, on the other hand, offers certain advantages over the classical isothermal method. First, this method eliminates the errors introduced by the thermal induction period, and second, it permits a rapid scan of the whole temperature range of interest. Several researchers have, therefore, employed the nonisothermal technique for oil shale pyrolysis (Campbell et al., 1978; Chen and Nuttall, 1979; Granoff and Nuttall, 1977; Haddadin and Mizyet, 1974; Herrell and Arnold, 1976; Rajeshwar, 1981). We decided to use both methods after Behnisch et al. (1980),who pointed out that a combined kinetic analysis of isothermal and nonisothermal TG measurements is an effective method for determining the most probable kinetic mechanism of polymer decomposition. Such a combined analysis was also performed by Campbell et al. (1978). In the present study, the weight loss data were analyzed by using three models. The kinetics of Moroccan oil shale are compared with those of Jordan and Colorado shale. 2. Experimental Section The oil shale used in this study originated from the Timhadit Site, Morocco. Chemical analysis of the oil shale is given in Table I. The samples were crushed and sieved to pass through a 200-mesh screen; the samples were used without further treatment. A Du Pont 951 thermogravimetric balance interfaced with a Du Pont 990 thermal analyzer was used to obtain the weight loss data under two sets of operating conditions: isothermal and nonisothermal. In the first set, the sample (25 mg) was introduced into the furnace maintained at a desired temperature. The isothermal oil shale decomposition was carried out at 325, 375, 410, 425,450, and 475 "C. In the second set of experiments, the sample (25 mg) was heated from ambient temperature to 900 "C at a linear heating rate in a dry stream of nitrogen. The heating rates of 1,2,5,10,20, and 50 OC/min were employed. The weight loss data and DTG curves were obtained in each case. The fraction of kerogen pyrolyzed, a, was defined by
where W, = initial weight of the sample (mg), W, = weight
of the sample at "t" minutes (mg), and W, = weight of the sample after complete pyrolysis of kerogen (mg). Each run was duplicated in order to minimize the error. 3. Kinetic Expression The general expression for the decomposition of a solid is given by Blazek (1973)
da/dt = kf(a)
(2)
where k = specific rate constant and f(a) = (1 - a ) for first-order reactions. If k is substituted in terms of activation energy and frequency factor, eq 2 can be rewritten as da/dt = Z exp(-E/RT)(l - a) (3) where Z = frequency factor (min-l), E = activation energy (J/g-mol), R = gas constant (J/(g-mol K)), and T = temperature (K). 3.1. Nonisothermal Study. For nonisothermal decomposition of solids, eq 3 can be modified by introducing heating rate as da dT = Z e x p -(1 - a ) (4) d T dt RET)
(
or d a / d T = (Z/b) exp(-E/RT)(l - a )
(5)
where b gives dTldt. In the present study, three models were used to evaluate the kinetic parameters: (1) Chen-Nuttall (1979), (2) Coats-Redfern (1964), and (3) Anthony-Howard (1976). Chen-Nuttall and Coats-Redfern models assume a single first-order rate equation to describe the decomposition reaction, while Anthony-Howard model assumes multiple parallel first-order reactions. Chen-Nuttall Model. In this method,. ea- 5 was rearranged and integrated to give 1 - a = exp{-(Z/b)L'exp(-E/RT)
or
I
1-a=exp--
dT)
Z RT2 bE+2RT
Rearranging eq 7 we get ZR b
1-a
E RT
Equation 8 is of the form y=mx+c where y = In
{
E+2RT T2
1 ln-)1 - a
m = -E/R c = In (ZR/b) x = 1/T
With the aid of repeated regression analysis, the values of slope and intercept can be computed. In the first step, an initial guess of E was made to calculate y (eq 10). Next, y was corrected for calculated values of E from eq 11, and m and c were recalculated from the regression. The iteration was continued until the desired accuracy of E and Z was satisfied.
Ind. Eng. Chem. Res., Vol. 26, No. 7 , 1987 1353 TEMPERATURE 473
673
I
I
OK
873
1073 I
f
0.1(nq/min)
i
. .
: : : :
I
(
, ,
I I
...-..................' 401
0
1
I
200
400
I
600
I
800
I
200
TEMPERATURE, *C
Figure 1. Nonisothermal TG curves for Moroccan oil shale a t various heating rates: (A) 1,(B) 2, (C)5, (D) 10, (E) 20, (F)50 K/min.
Coats-Redfern Model. Coats and Redfern (1964) developed a graphical method to determine the kinetic parameters for solid decompositions from nonisothermal data. The kerogen decomposition was considered to be first order in kerogen concentration as suggested in the literature (Rajeshwar, 1981; Campbell et al., 1978). According to this model, a plot of -[ln (-[ln ((1- a ) / T 2 ) ] ) ] against 1/T should result in a straight line of slope E / R . Anthony-Howard Model. Anthony and Howard (1976) proposed a model to explain the coal devolatization mechanisms. This model uses a Gaussian distribution function to account for the range of activation energies expected with such a mechanism. Recently, Nuttall et al. (1983) employed this model for oil shale decomposition in the form (for detailed derivation see Nuttall et al., 1983)
TEMPERATURE
A plot of -[ln (1- a ) ] against t yields a straight line with slope equal to k. The values of t used here include the induction period t o as suggested by Braun and Rothman (1975) and Pan et al. (1985). 4. Results and Discussion Figure 1presents the weight loss data for Moroccan oil
shale as a function of temperature at various heating rates. The important features of this figure are the following: (1) The total extractable kerogen content in Moroccan oil shale represents about 9 wt % of the total shale weight. The moisture content estimated from weight loss results below 200 "C corresponds to about 0.5-0.75%. The carbonate decomposition commenced at 525 "C or above depending upon the heating rate and gave a weight loss of about 25-26%. (2) Another observation concerns the effect of heating rate on total weight loss. A complete decomposition of kerogen was affected at 500 "C with a
€00
800
OC
Figure 2. DTG curves for Moroccan oil shale a t 2, 5, 20, and 50 K/min.
heating rate of 1 "C/min, while only 50% of total kerogen was found to decompose at 500 "C with a heating rate of 50 "C/min. In the latter case, the complete decomposition occurred at 600 "C. This difference is due to shorter exposure time to a particular temperature at faster heating rates. Figure 2 illustrates typical DTG curves corresponding to the data shown in Figure 1 for various heating rates. The salient features of this figure are as follows: (1)There is a shift in rate maximum toward higher temperature with the increase in heating rate from 5 to 50 "C/min in accordance with the literature data (Herrell and Arnold, 1976; Chen and Nuttall, 1979). (2) The presence of two clearly defined maxima in the temperature range 300-600 "C indicates that two reactions or processes occur consecutively in this temperature range. Actually, Allred (1966) has proposed a series mechanism for the pyrolysis of oil shale: kerogen
3.2. Isothermal Study. For isothermal decomposition of kerogen, eq 2 was integrated to give -[ln (1- a ) ] = k t (15)
400
kl
El,2,
gas bitumen carbon
k2
E2,Z,
gas oil coke
(16)
The present results are in accordance with his model and lead us to believe that the first peak (lower temperature) corresponds to the decomposition of kerogen to bitumen, while the second peak represents the decomposition of bitumen into oil and gas. These results are in agreement with those of Rajeshwar (1981). (3) The rate of weight loss increased by a factor of 3 with an increase in heating rate from 2 to 5 "C/min. However, a further increase in heating rate (5-50 "C/min) did not accelerate the decomposition rate appreciably. Figure 3 presents the isothermal weight loss data as a function of time in the temperature range 325-475 "C. During the initial period of these isothermal runs, less than 1 min was required to attain steady-state temperature. The use of small sample size (25 mg), small diameter furnace, and thin platinum foil boat aided in obtaining steady-state conditions in short time. During the temperature-induced time lag period, very small weight losses were recorded; at high temperatures (>450 "C), it amounted to 12-1470 of total weight loss. These results are very similar to those obtained by Pan et al. (1985).
1354 Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987 Table 11. Kinetic Parameters for the Nonisothermal and Isothermal DecomDosition of Moroccan Oil Shale activation energy, frequency factor, kJ/mol L/min ~~~~
method Coats-Redfern Chen-Nuttall isothermal
E, 40.2 43.9 38.4
E2 58.1 58.3 62.3
21
2 2
0.14 X lo3 0.11 X l o 3
0.19 X lo4 1.06 X IO4
5. Kinetic Analysis 5.1. Nonisothermal Study. In this study, two possible mechanisms were tested: (1) a single-step mechanism involving kerogen products at 300-500 "C, and (2) a two-step mechanism represented by eq 16. For the twostep mechanism, TGA data have been divided into two regions: region 1corresponding to the first step in a temperature range of 300-375 "C, and region 2 corresponding to the second step a t 375-500 "C. The coats-Redfern and Chen-Nuttall models were used to evaluate the kinetic parameters for both mechanisms, while the Anthony-Howard model was used to correlate pyrolysis data for the single-step mechanism at 300-500 "C. Two techniques, namely, iterative linear regression and the Lavenberg-Marquardt nonlinear least-squares regression, were employed to find the optimum values of E and Z in each of these equations. Two-step Mechanism. The Coats-Redfern analysis was performed by plotting -(ln [-On (1- a ) ) / T 2 ]vs. ) 1/T; the kerogen decomposition was assumed to be first order in kerogen concentration. Two sets of lines with different slopes (Figure 4) clearly indicate that two simultaneous reactions occur in the temperature range 300-600 "C. Activation energies were calculated from the slopes of these lines for each heating rate. The first reaction occurring at low temperature proceeds with an average activation energy of 40 kJ/mol, while the second reaction which starts at temperatures of 350 "C or above has an average activation energy of 58.1 kJ/mol. The fact that the variation in heating rate did not alter the slopes of these lines appreciably is in accordance with the findings of Rajeshwar (1981). The values of activation energy and frequency factor computed by using ChenNuttall and Coats-Redfern models are tabulated in Table 11. These models yield almost identical results. Single-Step Mechanism. The kinetic parameters and absolute average deviation (AAD) computed by using the Chen-Nuttall model/iterative linear regression technique, Coats-Redfern/iterative linear regression, and AnthonyHoward/Lavenberg-Marquardt nonlinear regression method are listed in Table 111. The kinetic parameters obtained by using Chen-Nuttall and Coats-Redfern models are almost identical. These values are, however, quite low as compared to those evaluated by Anthony-Howard equation. It can be noted from Table I11 that the values for activation energy and frequency factor as determined by the Anthony-Howard model are much more reasonable than those predicted by the other two models.
-
0
2
6
4
I 0 1 2 1 4
8
TIME (MINI
Figure 3. Percent weight loss for Moroccan oil shale at various temperatures (isothermal TG curves).
I 1.1
I
I
1.3
1.7
I ,I
IO'/T
Figure 4. Analysis of the TG data for Moroccan oil shale by the Coats-Redfern method. Heating rates: (0) 1, (A)2, (0) 5, ( 0 )10, (A)20, and (m) 50 K/min.
5.2. Isothermal Study. Figure 5 shows the plot of (1 as suggested by eq 15. Nearly straight lines, indicative of first-order kinetics, are obtained down
- a) vs. time plotted
Table 111. Kinetic Parameters and AAD Values for Single-Step Mechanism Chen-Nuttall Coats-Redfern heating E, 2, E, 2, rate, "C kJ/mol L/min AAD kJ/mol L/min 1 50.1 150.6 0.050 50.7 176.4 2 47.1 116.7 0.044 47.7 139.7 5 48.1 267.4 0.054 48.7 317.7 10 43.6 138.3 0.045 44.3 170.8 20 38.3 79.8 0.034 39.2 105.1 50 31.6 37.7 0.054 32.9 56.4
Anthony-Howard AAD 0.050
0.044 0.054 0.045 0.034 0.053
E, kJ/mol 90.1 90.4
Z,L/min 3.2
AAD 0.038
5.8
0.046
X
Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987 1355 Table IV. SDecific Rate Constants for Decomposition of Oil Shales (K/min-I) ~~
nonisothermal data for present study" a t 20 OC/min 50 OC/min
temp, OC 300 325 350 375 400 425 450 475 500 510
0.026 0.037 0.050 0.067 0.088 0.113 0.172 0.238 0.321
present study
0.050 0.065 0.085 0.107 0.134 0.170 0.200 0.236 0.400
0.05 0.09 0.15 0.20 0.31 0.41 0.63
isothermal data Jordon shalec Colorado shale (Haddadin and (Braun and Mizyet, 1974) Rothman, 1975)
0.050 0.070
0.02 0.06 0.18 0.48 2.90
0.113
0.033 0.053 0.218 0.390
(400)* (415) (430) (445)
0.143
o.05F 0
1985
0.020
k values calculated by using Coats-Redfern model. *Number in the parentheses gives the temperature in
0
Pan et al..
OC.
I
I
I
1.3
1.4
1.5
1.6
1
l000fT
Figure 6. Arrhenius plot in data. 2
.
4
6
0
8
IO3
12
13 4
T I M E (MIN)
Figure 5. Analysis of isothermal TG data (from Figure 3).
to a fractional weight change that depends upon temperature (Pan et al., 1985). Specific rate constants were evaluated from the slopes of the lines. Arrhenius plot (Figure 6) was used to determine activation energies and frequency factors. The kinetic parameters are compared with those obtained under nonisothermal conditions in Table 11. There is a close agreement between the values computed by eq 15 for isothermal runs and those obtained by using Coats-Redfern and Chen-Nuttall models for nonisothermal runs. The low activation energy (E,) for kerogen decomposition to bitumen compares favorably with the literature (Braun and Rothman, 1975; Rajeshwar, 1981) and indicates that the decomposition of kerogen to bitumen involves the breaking of relatively weak chemical bonds (Braun and Rothman, 1975). Haddadin and Mizyet (1974) and Haddadin and Tawarah (1980) also obtained low activation energy values and proposed a mechanism involving the diffusion of organic matter through the carbonate matrix as a rate-limiting step. The nonisothermal and isothermal data on specific rate constants are compared in Table IV. In the temperature range 300-425 "C, the isothermal data are comparable to nonisothermal data obtained at a heating rate of 50 "C/ min. This indicates that at lower temperatures, rates of isothermal decompositionare equivalent to those obtained under nonisothermal conditions with high heating rates (>50 OC/min). However, at temperatures exceeding 425 "C, k values for isothermal runs are 1.5 times higher than
In k vs. 1000/ T using isothermal TG
those obtained under nonisothermal runs. We feel that the most plausible explanation (among other explanations) for this discrepancy should be associated with the weight lost during induction period in the isothermal runs. At temperatures of 450 OC and above, about 12-14% weight loss was recorded during the thermal induction period. Although this initial period was less than 1 min, 12-14% weight was lost (under nonisothermal conditions) probably at a heating rate exceeding 450 OC/min. This may be responsible for the higher 12 values in isothermal runs. Table IV also compares the isothermal kintics ( k values) for Moroccan shale pyrolysis with the literature values for the decomposition of Jodan (Hadaddin and Mizyet, 1974) and Colorado shale (Braun and Rothman, 1975; Pan et al., 1985). This comparison shows that at temperatures lower than 425 OC, the decomposition rate for Moroccan shale is 3 times faster than that for Jordan or Colorado shale. At 425-475 "C, the rate is about 4 times faster than that for Jordan shale and almost equvalent to that for Colorado shale.
Conclusions The combined use of nonisothermal and isothermal TG measurements has shown that the thermal decomposition of Moroccan oil shale involves two consecutive reactions with bitumen as an intermediate. Both reactions follow first-order kinetics. Among the three models used in this study, AnthonyHoward model yields lower deviation and, thus, provides a better fit of the data. The decomposition rate of Moroccan oil shale is higher than Jordan or Colorado shale at temperatures lower than 425 "C. At 425-475 "C, it remains higher than that for
Ind. Eng. Chem. Res. 1987,26, 1356-1362
1356
Jordan shale but is equivalent to that for Colorado shale. A t temperatures higher than 475 "C, pyrolysis rate of Moroccan shale is lower than Colorado shale. Literature Cited Allred, V. D. Chem. Eng. Prog. Symp. Ser. 1966, 62(8), 55. Anthony, D. B.; Howard, J. B. AIChE J . 1976,22, 625. Behnisch, J.; Schaff, E.; Zimmermann, H. Thermochim. Acta 1980, 42, 65. Blazek, A. Thermal Analysis; Tyson, J. F., Transl. Ed.; Van Nostrand-Reinhold: London, 1973. Braun, R. L.; Rothman, A. J. Fuel 1975,57, 129. Campbell, J. H.; Koskinas, G.; Stout, N. Fuel 1978, 57, 372. Chen, W. J.; Nuttall, H. E. Paper presented at the 86th AIChE National Meeting, Houston, TX, 1979.
Coats, A. W.; Redfern, J. P. Nature (London) 1964, 201, 68. Dericco, L.; Barrick, P. L. Ind. Eng. Chem. 1956, 48, 1316. Galan, M. A,; Smith, J. M. AIChE. J . 1983, 29, 604. Granoff, B.; Nuttall, H. E. Fuel 1977, 56, 234. Haddadin, R. A.; Mizyet, F. A. Ind. Eng. Chem. Process Des. Deu. 1974, 13(4), 332. Haddadin, R. A.; Tawarah, K. M. Fuel 1980, 59, 539. Herrell, A. Y.; Arnold, C., Jr. Thermochim. Acta 1976, 17, 165. Hubbard, A. B.; Robinson, W. E. Rep. Invest.-L1.S. Bur. Mines 1950, 4744, 1. Nuttall, H. E.; Guo, T.-M.; Schrader, S.; Thakur, D. S. ACS Symp. Ser. 1983, 230, 269. Pan, Z.; Feng, H. Y.; Smith, J. M. AZChE. J . 1985, 31, 721. Rajeshwar, K. Thermochim. Acta 1981, 45, 253. Received for review March 10, 1986 Accepted March 16, 1987
Stability of Tubular and Autothermal Packed Bed Reactors Using Phase Plane Analysis Richard W. Chylla, Jr., Raymond A. Adomaitis, and Ali Char* Department of Chemical Engineering, Illinois Institute of Technology, Chicago, Illinois 60616
The regions of stability and parametric sensitivity of countercurrent reactor/heat exchangers are determined explicitly in the plane of inlet feed temperature-inlet coolant temperature. The concept of phase plane analysis is generalized to include all orders of reaction rate expressions, a broader range of system parameters, and is extended to the case of autothermal reactors. An industrial hydrocarbon oxidation reactor model and an autothermal CO oxidation reactor model have been used to illustrate and to evaluate the analysis method. The approach presented here is appealing since the region of safe inlet temperatures is determined explicitly and the region of safe operation can be optimized with respect to the reactor design parameters. 1. Introduction
Industrial packed bed reactors serve as the workhorse of the chemical industry, and yet there is much progress to be made in terms of obtaining reliable models for reactor design and control purposes. Autothermal reactors in particular suffer from large parametric sensitivities due to the inherent feedback mechanism. A large step fonvard in the design and operation of these reactors would be a diagnostic method to determine the region of safe operating temperatures for both the inlet coolant and the reactants. Accurate models of reactor temperature behavior are often complex and even numerical solution on a digital computer can be unwieldy and time consuming. Perhaps more importantly, the effect of a particular design parameter on the stability of the reactor is not always obvious and may often require repeated solutions. A major improvement in any control effort is to maximize the open-loopregion of stability. This work is aimed at providing a tool for the engineer designing the reactor and/or the control system to determine a priori how stable a given reactor configuration is and how the design parameters will affect the size of this region of stability. There have been numerous studies in the literature which offer criteria with which one may determine whether or not a set of inlet conditions leads to stable operation (van Welsenaere and Froment, 1970; Oroskar and Stern, 1979; Froment and Bischoff, 1979; Bilous and Amundson, 1956). The concept of using phase plane analysis to determine regions of safe operation has been proposed by Akella and Lee (1983) for reactors cooled by a countercurrent fluid, especially liquids. In this context, phase
* Author to whom correspondence should be addressed. 0888-5885/87/2626-1356$01.50/0
plane analysis is used to denote the reactant temperaturecoolant temperature plane. Their algorithm is limited to positive-order reactions, and autothermal operation was not considered. The phase plane analysis in this paper: (1) extends analysis of Akella and Lee (1983) to a more general class of reactions, (2) permits the designer of an autothermal reactor to know how safe a set of inlet temperatures are for a given set of operation conditions, and (3) determines the necessity of a feed-effluent heat exchanger for autothermal operation. 2. Phase Plane Analysis of T u b u l a r Reactors The vast majority of industrial reactions take pldce in fixed bed catalytic reactors, and many of the bulk &emicals such as ammonia, methanol, and functionalized aromatics are produced by exothermic reactions. This requires large amounts of heat to be removed to prevent thermal runaway, undesirable side reactions, or unfavorable equilibrium conditions. In many reactors, this heat is removed by an external medium such as water or heat stable molten salts. In large scale operations, it is economical to use the heat of reaction to preheat the feed to desired reaction temperature. Since the coolant eventually enters the reactor as the feed stream, these reactors display high levels of parametric sensitivity and very narrow windows of safe operation. The following discussion reviews the development of an algorithm to produce a phase plane for a general countercurrent reactor/heat-exchanger configuration. Autothermal reactors are then treated as a subset of these reactors. Industrial exothermic reactions are commonly carried out in tubular packed beds filled with catalyst and cooled by flowing a cooling medium countercurrently. The 0 1987 American Chemical Society
