Ind. Eng. Chem. Res. 2004, 43, 3723-3730
3723
Time to Runaway in a Continuous Stirred Tank Reactor Dale D. Slaback and James B. Riggs* Department of Chemical Engineering, Texas Tech University, Lubbock, Texas 79409-3121
Thermal runaway is a significant safety concern in chemical plants with exothermic reactors. Understanding the time scale on which thermal runaway occurs for these systems can be an important factor for their safe and efficient operation. On the basis of the knowledge of time to runaway, procedures can be implemented to reduce the probability of occurrence of a thermal runaway. This paper develops correlations for time to runaway for a continuous stirred tank reactor (CSTR). The system in this study involves a reactor with a single, irreversible, firstorder, exothermic reaction. Correlations for time to runaway are developed for both adiabatic and nonadiabatic reactors. These correlations are based on results obtained by simulating the reactor system for a range of values of the dimensionless parameters that define the system. The system time constant, developed using an eigenvalue analysis, is included in the time to runaway correlations. In addition, a new definition of runaway, based on the dynamic temperature profile, is developed. I. Introduction Exothermic reactors are widely used in the chemical process industries, and thermal runaway is a significant safety concern for these exothermic reactors. Thermal runaway results because an increase in the reactor temperature causes an increase in the heat generated by the exothermic reaction. If the heat removal equipment is not adequate, sharp increases in the reactor temperature can result, triggering an explosion for certain cases. Some common reactors that are prone to thermal runaway are ethylene oxide reactors, polymer reactors, and hydrocarbon nitration reactors. Because of the risk of thermal runaway, many industrial reactors are operated in an overly conservative, suboptimal manner. As a result, understanding thermal runaway is of significant industrial importance. Thermal explosion theory was first applied to chemical reactor systems by Semenov.1 Bilous and Amundson2 furthered the study in this area, introducing the concept of thermal runaway. Morbidelli and Varma3-5 and van Welsenaere and Froment6 studied the parametric sensitivity of tubular and fixed-bed reactors with a single reaction. McGreavy and Adderley7 added the effects of particle mass. Reactors with multiple reactions were studied by Hosten and Froment,8 Henning and Perez,9 and Morbidelli and Varma.10 The effect of coolant temperature was investigated by Soria Lopez et al.,11 Hosten and Froment,8 Bauman et al.,12 and Henning and Perez.9 Emig et al.13 verified the accuracy of several of these criteria experimentally. Kimura and Levenspiel14 developed a graphical procedure to determine the regions of parametric sensitivity. Chemburkar et al.15 expanded the parametric sensitivity studies to continuous stirred tank reactors. The whole of this research has focused on developing steady-state criteria for reactor runaway. Specifically, these studies have identified regions in which small changes in the reactor operating conditions lead to large steady-state changes in the reactor temperature. Because this research has been performed on steady-state * To whom correspondence should be addressed. Tel.: (806) 742-1765. Fax: (806) 742-1765. E-mail:
[email protected].
systems, the results indicate only when the system is prone to thermal runaway. Once the system parameters reach values indicating runaway, however, the parametric sensitivity results provide no indication of how long it takes for the runaway process to occur, the time to runaway. In this study, time to runaway is calculated from the time at which the system parameters reach values indicating that runaway will occur. In many industrial applications, however, thermal runaway cannot be detected until a significant amount of time after it becomes imminent. For this reason, the exact time to runaway cannot be predicted. However, having an estimate for the time to runaway for a given system can provide an indication of the severity of the problem associated with thermal runaway. For systems with long times to runaway, dealing with a thermal runaway is relatively simple. The system can be operating near the sensitive region, and if the conditions in the reactor do vary such that runaway becomes imminent, the long time to runaway allows the system to be corrected before the runaway occurs. For systems with short times to runaway, however, dealing with thermal runaway is considerably more difficult. For these systems, the reactor must be operated away from the sensitive region to prevent a fast, sometimes nearly immediate, runaway, and specific safety overrides must be implemented. The focus of this study is to develop an estimate of the time to runaway, giving an indication of the severity of the problem associated with thermal runaway. The remainder of this paper considers the modeling of the time to runaway for a CSTR with a first-order reaction. II. Reactor Modeling and Simulation A dynamic model is used in this study to calculate the time to runaway for a wide variety of operating conditions. The system considered in this study is a continuous stirred tank reactor in which a single, irreversible, first-order, exothermic reaction occurs. The modeling of this system neglects the latent heat and radiation effects. With these assumptions, the material and energy balance equations for this system are
10.1021/ie0306574 CCC: $27.50 © 2004 American Chemical Society Published on Web 01/15/2004
3724 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
dCA ) -r dt Fcp
Table 1. Parameter Ranges for Development Data
(1)
dT v UA ) Fcp(Tinlet - T) + r(-∆Hrxn) (T - Tc) dt V V (2) CA(0) ) CA0,
T(0) ) T0
where
r ) k0CAe-E/RT The material and energy balance equations were converted into dimensionless form by defining reference values for concentration, temperature, and time. The reference values chosen were the inlet concentration and temperature and the reactor volume-to-flow rate ratio. In addition to the reference values, several dimensionless groups were defined for use in the material and energy balances: the Damko¨hler number, dimensionless activation energy, dimensionless heat of reaction, and Stanton number. Also, the dimensionless coolant temperature and dimensionless reactor inlet temperature were defined. The definitions of these groups are as follows
k(Tinlet)V v
(3)
E RTinlet
(4)
(-∆Hrxn)CA,inlet FcpTinlet
(5)
St )
UA Fcpv
(6)
Θc )
Tc Tinlet
(7)
Da ) γ)
B)
Θinlet )
Tinlet )1 Tinlet
(8)
Substituting these dimensionless groups into the material and energy balance equations allows them to be expressed more succinctly as
dXA ) Da(1 - XA)eγ(Θ - 1)/Θ dτ
(9)
dΘ ) (1 - Θ) + BDa(1 - XA)eγ(Θ - 1)/Θ - St(Θ - Θc) dτ (10) XA(0) ) 0,
Θ(0) ) Θ0
The material and energy balance equations can be further simplified by assuming that the reactant is present in large excess, indicating that the conversion is equal to zero. If the reactor conversion were not assumed to be zero, the reactor temperature would increase more slowly, leading to longer times to runaway. Therefore, the times to runaway predicted in this study are minimum values for each system. For systems
dimensionless parameter
minimum value
maximum value
heat of reaction × Damko¨hler number activation energy initial reactor temperature Stanton number coolant temperature
0.05 15 1.1 0 0.9
0.25 19 1.3 20 1.1
operating with a large excess of reactant, the actual time to runaway will be close to this minimum value. For systems without a large excess of reactant, however, the runaway will take significantly longer than the predicted minimum time to runaway. The system can now be defined entirely by the energy balance, noting that the conversion term can be eliminated. This simplified energy balance equation is used to represent the system for the remainder of this study. The simplified energy balance contains five parameters, defined in eqs 3-7, that define the behavior of the system. The dimensionless heat of reaction and the Damko¨hler number appear in the energy balance equation only as the product of the two. For this reason, the product of the dimensionless heat of reaction and the Damko¨hler number was treated as a single parameter. With this convention, the energy balance equation contains four parameters (product of dimensionless heat of reaction and Damko¨hler number, dimensionless activation energy, Stanton number, and dimensionless coolant temperature) that define the behavior of the system. In addition to these four parameters, the initial dimensionless temperature of the reactor also affects the system behavior. A set of development data was generated by integrating the dimensionless energy balance equation over time using the LSODE16 integrator. These development data were used to develop correlations for time to runaway. The values of the five dimensionless parameters were varied over a range consistent with the data provided by van Welsenaere and Froment6 and Soria Lopez et al.11 The ranges for each of the dimensionless parameters are summarized in Table 1. The maximum value of the initial reactor temperature was chosen to prevent a significant number of nearly immediate runaways. Similarly, the minimum value of the coolant temperature was chosen to prevent an excessive number of stable reactor results. The maximum value of the coolant temperature coincides with the minimum initial reactor temperature so that no combination of parameters gives a reactor that is heated by the coolant. Within these ranges, each parameter was assigned five possible values: the maximum, the minimum, and three equally spaced values between these extremes. The development data set includes the results from simulations using every combination of the possible values for each of the system parameters. As a result, the development data set contained 3125 sets of simulation results, with 1952 of these sets representing runaway reactors. The parameter ranges used were such that some combinations resulted in a stable reactor, others yielded long times to runaway, and still others led to runaway in very short times. Additionally, the special case of an adiabatic reactor, in which the Stanton number is equal to zero, was considered. In general, the temperature profiles produced by the simulations that result in thermal runaway exhibit a hyperbolic functionality. The profile begins with a relatively flat region. After some time, the profile begins to curve sharply upward. Shortly
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3725
Figure 1. Typical dimensionless temperature profiles. Table 2. Parameters for Typical Simulation Profiles profile dimensionless parameter heat of reaction × Damko¨hler number activation energy initial reactor temperature Stanton number coolant temperature
1
2
3
4
5
6
0.2 0.2 0.2
0.25 0.25 0.15
15 15 18 1.3 1.2 1.15 0 0 10 1.0
15 1.15 10 1.0
17 1.1 10 1.0
16 1.15 15 0.95
after the curvature begins, the temperature profile becomes nearly vertical. Although the time scale and sharpness of the curve vary from run to run, the overall behavior of all runs is similar. Figure 1 shows some typical dimensionless temperature profiles for parameter sets resulting in both short and longer times to runaway. This figure also indicates the onset of thermal runaway for each profile. As will be defined later, the onset of thermal runaway corresponds to a slope of unity in the dimensionless temperature profile. The parameters for each of these profiles are listed in Table 2. Because the first two profiles are for adiabatic simulations, the coolant temperature is not used in these simulations. Although these different sets of parameter values result in different times to runaway, the functional behaviors of all of the simulations are similar. Inspection of the dimensionless temperature profiles shown in Figure 1 leads to some conclusions about the time to runaway. Profiles 1 and 2 are for systems that differ only in the initial dimensionless reactor temperature. As can be seen in the figure, thermal runaway occurs at the same dimensionless temperature for each of these two systems. This indicates that, although the time to runaway is dependent on the initial dimensionless temperature, the dimensionless temperature at which the runaway occurs is independent of the initial dimensionless temperature. Also, profiles 3, 4, and 6 all have the same initial dimensionless temperature. The times to runaway and the dimensionless temperatures at which runaway occurs are significantly different among these profiles, indicating that the additional operating parameters significantly affect both the time to runaway and the dimensionless temperature at which runway occurs. Finally, profile 5 begins at a lower initial dimensionless temperature than profile 6. However, thermal runaway occurs sooner in the system represented by profile 5. This comparison indicates that the time to runaway has a stronger dependence on the additional operating conditions than on the initial dimensionless temperature. III. Definition of Time to Runaway The runaway condition has been defined in various ways in the classical parametric sensitivity literature.
Chambre´17 and Barkelew18 introduced the use of isoclines in the determination of the runaway criterion. This approach provides a graphical method for determining runaway conditions in a reactor. Several other studies have used modifications of this approach, defining the runaway condition by locating maxima in some parameter of the reaction system. All of this work, however, was performed using a steady-state analysis of the system. To define runaway for a dynamic system, a particular characteristic of the dimensionless temperature profile was used. Because this characteristic must exist for all runaway simulations, the first candidate for defining the onset of thermal runaway is to select a specific value of the dimensionless temperature. With the assumptions made in the process model, two possible simulation results exist: stable and runaway. For stable simulations, a maximum achievable dimensionless temperature exists. This maximum temperature could possibly provide an indication of the onset of thermal runaway. Although Semenov1 reports this maximum to be unity for all reactive systems, some equations relating this maximum to the reaction order can be developed from steadystate runaway criteria developed in the parametric sensitivity studies. Using the Adler and Enig19 criterion, the maximum achievable dimensionless temperature in a stable exothermic reactor is given by
ΘM ) 1 + xn
(11)
where n is the overall order of the reaction. Similarly, using the Thomas and Bowes20 criterion, the maximum achievable dimensionless temperature can be expressed as
ΘM )
2 - n + xn2 + 4n 2
(12)
For the first-order system studied in this work, these expressions give maximum achievable dimensionless temperatures of 2 and 1.618, respectively. The dimensionless temperature profiles shown in Figure 1 indicate that the system begins to run away well before the dimensionless temperature reaches either of these values. Upon further inspection of the simulation results, no single dimensionless temperature seems to indicate the onset of runaway for all of the parameter sets. Therefore, another definition of runaway must be developed. The slope of the dimensionless temperature profile can also be used to determine the onset of thermal runaway. Each dimensionless temperature profile exhibits a slow rise in dimensionless temperature followed by sudden curvature and a rapid increase of the dimensionless temperature. Because of this rapid increase, the dimensionless temperature becomes quite large very shortly after the curvature begins. The onset of thermal runaway is defined in this study as the point at which the curvature in the dimensionless temperature profile is the largest. For all of the dimensionless temperature profiles, this point occurs when the dimensionless temperature profile reaches a 45° angle with respect to the dimensionless time axis. The slope of the dimensionless temperature profile at this point is equal to 1. Therefore, the onset of thermal runaway was defined in this study as the temperature at which the slope of the dimensionless temperature profile equals 1.
3726 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
To determine the dimensionless temperature at which thermal runaway occurs, the critical value for the slope of the dimensionless temperature profile was substituted into the dimensionless energy balance to give
1 ) (1 - Θ) + BDaeγ(Θ - 1)/Θ - St(Θ - Θc) (13) For each set of system parameters, eq 13 is a single equation in a single variable, the dimensionless temperature. Therefore, solving this equation gives the dimensionless temperature at which thermal runaway begins. Because the equation is nonlinear, multiple solutions exist, and the solution technique and initial estimate of the temperature can affect the result. In this study, the dimensionless temperature at the onset of runaway was calculated using Newton’s method.21 The maximum achievable dimensionless temperature as predicted from the Thomas and Bowes criterion was used as the initial estimate. In general, the runaway temperature can be expressed as a function of the system parameters such that
Θra ) f(γ,B,Da,St,Θc)
(14)
Thermal runaway can now be identified for each system by calculating the appropriate root and monitoring the system temperature. This approach produces a runaway temperature reliant upon the entire set of system parameters, not just the reaction order. The time to runaway can then be determined as the amount of time for which the reactor operates until this temperature is achieved. IV. Analytical Calculation of Time to Runaway Ideally, the dimensionless time to runaway would be determined directly from the dimensionless energy balance equation. Integration of the energy balance equation would give an explicit relationship between dimensionless time and temperature. To integrate the energy balance equation, separation of variables must be performed. Because the right-hand side of the dimensionless energy balance does not contain the dimensionless time, separation of variables is straightforward. Because of the exponential term in the denominator of the separated equation, an analytical form of the integral does not exist. The exponential, therefore, must be approximated. Because the exponential is defined by an infinite series, a truncation of this series is used as the exponential approximation. The simplest approximation that can be used is the first-order approximation. Substituting this approximation into the separated energy balance equation allows the equation to be integrated to give a closed-form expression for time to runaway. This function, however, does not give accurate values for the time to runaway. The accuracy of the analytical solution could be improved by using more terms in the exponential approximation. Unfortunately, a more accurate exponential approximation results in an energy balance equation that does not lend itself to development of a closed-form expression for the time to runaway. Using the runaway temperature in the exponential approximation could also increase the accuracy of the solution. However, the improvement gained by using this temperature is minimal.
V. Empirical Modeling of Time to Runaways Adiabatic Case Because an accurate analytical solution for time to runaway could not be generated, a simple yet accurate numerical approximation was sought. This approximation was initially developed for the special case of the adiabatic reactor. As the Stanton number is equal to zero for this case, the system is now described by only three parameters: the product of the dimensionless heat of reaction and Damko¨hler number, the dimensionless activation energy, and the initial dimensionless reactor temperature. First, the development of a simple linear regression relating the time to runaway and these three parameters was attempted. Because the time to runaway varies inversely with each parameter, the regression correlates the time to runaway with the reciprocals of the system parameters. Unfortunately, no acceptable regression could be produced to relate the time to runaway to the system parameters alone. Because the simple regression did not give acceptable results, a regression with additional information must be developed. The information included in this regression must provide some insight into the time frame on which the system operates. One such characteristic of the system is the time constant as calculated from an eigenvalue analysis of the system. Although the eigenvalue is based on a linearization of the nonlinear equation, its inclusion in the regression might produce acceptable results. The eigenvalue of the system is calculated by differentiating the right-hand side of the energy balance equation at the initial conditions, giving
1 λ ) -1 + BDaγ 2eγ(Θ0 - 1)/Θ0 - St Θ0
(15)
The time constant for the system is equal to the reciprocal of the eigenvalue. Inverting the eigenvalue and rearranging results in the expression
τλ )
Θ02 BDaγeγ(Θ0 - 1)/Θ0 - Θ02(1 + St)
(16)
Next, a regression including the time constant was developed using the data set previously described in Table 1. This regression was developed using the leastsquares approach employed in the Microsoft Excel data analysis package. Because time to runaway varies inversely with each of the system parameters, the reciprocals of the parameters were also included in the regression. The results of this regression indicated that two of the parameters (dimensionless heat of reaction and initial dimensionless reactor temperature) affect the time to runaway only through the system time constant. Therefore, these parameters were removed from the regression. The time to runaway data was again regressed, but this time as a function of the time constant and dimensionless activation energy only. The resulting time to runaway for the adiabatic case can be expressed as
1 τA ) 1.12τλ - 0.730 + 0.006 30 γ
(17)
The goodness of fit of this regression was assessed using several statistics. Judging from these statistics,
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3727 Table 3. Statistical Parameters Used to Determine Goodness of Fit regression
multiple correlation coefficient
mean squared error
slope standard error
intercept standard error
adiabaticsdevelopment adiabaticstesting nonadiabatic, standardsdevelopment nonadiabatic, near criticalsdevelopment nonadiabatic, similar temperaturessdevelopment nonadiabaticstesting
0.9991 0.9975 0.9684 0.9673 0.9858 0.9538
3.00 × 10-5 4.94 × 10-4 8.04 × 10-5 5.87 × 10-4 7.46 × 10-5 3.50 × 10-4
0.0017 0.0086 0.052 0.0642 0.0281 0.046
0.0002 0.001 0.0015 0.0063 0.0018 0.005
Table 4. Parameter Ranges for Test Data dimensionless parameter
minimum value
maximum value
heat of reaction × Damko¨hler number activation energy initial reactor temperature Stanton number coolant temperature
0.03 14 1.0 0 0.9
0.3 25 1.4 25 1.1
systems are less susceptible to severe consequences than the systems with short times to runaway. Therefore, operating a reactor in a manner that increases the time to runaway would be important to loss prevention. Judging from the form of the correlation developed, the time to runaway can be increased in several ways. First, increasing the reactor inlet temperature while maintaining the same reactor temperature leads to an increase in the time to runaway. This could be achieved by manipulating the duty of a feed heater or cooler. Conversely, the reactor temperature can be decreased while maintaining a constant inlet temperature. Judging from the forms of certain dimensionless parameters, changes in the reacting fluid flow rate can lead to either increases or decreases in the time to runaway. However, at the normal operating conditions, increasing the reacting fluid flow rate results in a longer time to runaway. VI. Empirical Modeling of Time to Runaways Nonadiabatic Case
Figure 2. Predicted vs simulated time to runaway in adiabatic reactors.
the correlation fits the data quite well, with a multiple correlation coefficient of 0.9991 and a mean squared error of 3.00 × 10-5. The values of the statistical goodness of fit parameters for each regression and data set in this study are summarized in Table 3. In addition to the data for which the correlation was developed, some additional parameter sets were used to verify the accuracy of the correlation. These parameter sets include values in the same range as the initial set as well as some values outside of this range. These data were generated by integrating the dimensionless energy balance equation over time for the parameter ranges given in Table 4. For this set of data, each parameter had three to five possible levels, and every parameter combination was simulated. As a result, 775 simulation results were obtained. Of these results, 616 sets represented runaway reactors and were included in the test data set. In general, the correlation predicts accurate values for these data as well. The exception to the proper fit is for initial temperatures below the range for which the correlation was developed. This region represents a system in which the initial reactor temperature is equal to the feed temperature, a situation that is not likely in practice. For initial temperatures above the range used in development, as well as for dimensionless activation energies and dimensionless heats of reaction (times the Damko¨hler number) both above and below the range used in development, the correlation fits well. Figure 2 shows the predicted and simulated values of the time to runaway of the adiabatic test data. Because systems with longer times to runaway can be corrected before the thermal runaway occurs, these
A regression analysis for the nonadiabatic case was also performed. As in the adiabatic case, the time constant obtained from the eigenvalue analysis of the system and all of the system parameters were included in this regression. The result of this regression is a correlation similar to that developed for the adiabatic case. However, for nonadiabatic systems, several predictions were found to have significant error, predictions that were larger or smaller than the simulated value. Each parameter set that produced a prediction of the time to runaway that was less than the simulation result contains the lowest value of the dimensionless heat of reaction parameter that results in runaway. That is, when the dimensionless heat of reaction parameter is decreased to the next lower value, the simulation results in stable operation. Consideration of how a heat of reaction parameter near the stability limit affects the dimensionless temperature profile validates this finding. For such a simulation, the initial slope of the temperature profile is nearly zero. Therefore, the temperature increases very slowly at the beginning of the simulation. Only after a significant incubation period does the exponential nature of the system begin to have a significant effect on the temperature profile. This phenomenon causes the simulated time to runaway to be longer than for cases not near this stability limit. Because this effect was not considered in the correlation, the correlation predicted a time to runaway less than the simulated value. The stability limit occurs when the derivative of the dimensionless temperature with respect to dimensionless time equals zero. The energy balance for this system at the stability limit is
0 ) (1 - Θ) + BDaeγ(Θ - 1)/Θ - St(Θ - Θc) (18)
3728 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
The energy balance equation is solved for the dimensionless heat of reaction (times the the Damko¨hler number) to give the critical dimensionless heat of reaction parameter. Because the critical value at the initial condition is needed, the dimensionless reactor temperature is taken at its initial value. The critical dimensionless heat of reaction parameter is calculated as
BDacritical )
(St + 1)Θ0 - Θin - StΘc eγ(Θ0 - 1)/Θ0
(19)
To use this critical value, the concept of being near the critical condition must be defined. In this study, values of the dimensionless heat of reaction parameter within 25% of the critical value were considered near critical. An analysis of the parameter sets resulting in time to runaway predictions greater than the simulation result showed that such values result when the initial reactor temperature and coolant temperature are close to one another. An analysis similar to that applied for near-critical conditions served to validate this result as well. Because the difference in the reactor and coolant temperatures is small, the reactor initially exhibits limited heat transfer. Only after some time, when the reactor temperature has increased, is the total heat transfer achieved. Because the heat transfer is initially low, less heat is removed, and the temperature increases more quickly than predicted by the time constant. As a result, the correlation predicts a value for the time to runaway that is higher than the simulation result. Analogously to the near-critical case, similar temperatures must be defined. Similar-temperature is defined in this study as a condition for which the dimensionless initial reactor temperature and dimensionless coolant temperature have an absolute difference of no more than 0.05. With the causes of the error identified, the nonadiabatic data were separated into three categories: standard, near critical, and similar temperatures. Each of these categories was considered separately in the analysis. The first category of data, standard, was regressed versus the input parameters as previously attempted. This reduced data set produced a correlation with a much better fit of the data. Again, some of the parameters were shown to affect the time to runaway only through the system time constant. The correlation obtained for the standard nonadiabatic systems is
1 1 τNA ) 1.14τλ - 0.698 + 0.0954 - 0.0773 (20) γ Θc As would be expected, this correlation does not fit the nonadiabatic data quite as well as the adiabatic regression fits the adiabatic data. However, the correlation still fits the data fairly well, with a multiple correlation coefficient of 0.9684 and a mean squared error of 8.04 × 10-5. The statistical goodness of fit parameters are summarized in Table 3. Once the correlation for the standard data was completed, the exceptions were considered. Rather than develop a new correlation for these data sets, correction factors for the standard correlation were developed such that the times to runaway for near-critical and similartemperature parameter sets are calculated as
Figure 3. Predicted vs simulated time to runaway for nonadiabatic reactors.
τNC ) τNA + CFNC
(21)
τST ) τNA + CFST
(22)
These correction factors were determined by regression of the errors resulting from the standard nonadiabatic correlation. In addition to the original parameter set, the regression for the near-critical correction factor also included the critical value of the heat of reaction parameter and the ratio of the actual parameter value to the critical value. The resulting near-critical correction factor is
1 1 1 CFNC ) -0.197τλ + 3.00 - 0.0114 + 0.707 γ Θc Θ0 BDa 0.419BDacritical - 0.384 - 0.198 (23) BDacritical In a similar manner, the regression for the similartemperature correction factor included the temperature difference as a parameter. The resulting similar-temperature correction factor is
1 1 1 CFST ) -0.811τλ + 2.42 + 4.22 - 3.70 + γ Θc Θ0 1 0.009 89 - 3.15(Θ0 - Θc) - 0.617 (24) BDa The near-critical and similar-temperature regressions also fit the data quite well, as indicated by the statistical goodness of fit parameters summarized in Table 3. The multiple correlation coefficients for the near-critical and similar-temperature correlations are 0.9673 and 0.9858, respectively. The mean squared errors for these correlations are 5.87 × 10-4 and 7.46 × 10-5, respectively. With the correction factors, the time to runaway for nonadiabatic reactors can be predicted for the test data set generated using the parameter values in Table 4. For the nonadiabatic test data, 83% of the simulation results were for standard conditions, 7% for near-critical operation, and 11% for similar-temperature operation. Figure 3 shows the predicted versus simulated times to runaway for the test data for standard, near-critical, and similar-temperature nonadiabatic reactors. Again, operating these reactors in a manner that increases the time to runaway is important in loss prevention. As in the adiabatic system, changes to the reacting fluid stream can increase the time to runaway. The coolant stream for a nonadiabatic system can also be manipulated to increase the time to runaway.
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3729 Table 5. Case Study Data6,11 parameter
value
units
heat of reaction activation energy inlet temperature reaction rate at inlet temperature volumetric flow rate reactor volume inlet reactant concentration reacting fluid density reacting fluid heat capacity initial reactor temperature overall heat-transfer coefficient heat-transfer surface area coolant temperature
-307 000 19.712 620 1.70 × 10-3 1.767 4.418 0.02 1.293 0.25 744 82.7 0.0345 620
kcal/kmol kcal/kmol K 1/h m3/h m3 kmol/m3 kg/m3 kcal/kg‚K K kcal/m2‚h‚K m2 K
Table 7. Average Operating Data for an Ethylene Oxide Reactor22 parameter
value
units
heat of reaction activation energy inlet temperature reaction rate at inlet temperature volumetric flow rate reactor volume inlet reactant concentration reacting fluid density reacting fluid heat capacity initial reactor temperature overall heat-transfer coefficient heat-transfer surface area coolant temperature
-230 659 137 500 3.49 × 10-4 1.08 × 10-4 0.01 0.5 13 1858.4 550 20 0.5 450
J/kmol kJ/kmol K 1/s m3/s m3 mol/L kg/m3 J/kg‚K K J/m2‚s‚K m2 K
Table 6. Case Study Dimensionless Parameters dimensionless parameter
value
heat of reaction × Damko¨hler number activation energy initial reactor temperature Stanton number coolant temperature
0.1299 16 1.2 4.9992 1
Increasing the total heat transfer leads to an increase in the time to runaway. Either decreasing the coolant temperature or increasing the coolant flow rate can achieve this condition. VII. Case Studies To demonstrate the predictive capabilities of the correlations developed in this study, two case studies were performed. The first case study was performed on the system described in the parametric sensitivity literature.6,11 The operating conditions for this system are reported in Table 5. Equations 3-7 were used to convert the raw data for the system into dimensionless form. The dimensionless parameters for the system considered in this case study are listed in Table 6. This system was simulated as before, and the time to runaway was determined from the runaway definition previously developed. The dimensionless time to runaway for this system was found to be 0.0587 (352 s). Next, the appropriate correlation was applied to calculate the time to runaway. Because the system is nonadiabatic, it must be classified as standard, similar-temperature, or near-critical operation. The initial dimensionless reactor temperature and dimensionless coolant temperature are not similar. In addition, calculation of the critical heat of reaction parameter indicates that the system is not near the critical condition. Therefore, this system represents a standard nonadiabatic reactor. Accordingly, the time to runaway for this reactor was calculated using eq 20. This equation gave a dimensionless time to runaway of 0.0516 (310 s) for the system. The error of this prediction is 0.0071 (42 s), corresponding to a 12% error. Despite the error in the calculated time to runaway, thermal runaway is correctly classified as a fairly difficult problem for this system. The time frame of the runaway is such that the process cannot merely be corrected once the runaway is identified. The process must be operated away from the sensitive regions and must be equipped with fast-acting overrides to combat thermal runaway if it does become imminent. The predictive capability of these correlations was also tested for reactions occurring in tubular reactors. Because the operating parameters for a tubular reactor
Table 8. Dimensionless Parameters for Ethylene Oxide Reactor dimensionless parameter
value
heat of reaction × Damko¨hler number activation energy initial reactor temperature Stanton number coolant temperature
0.3099 32.956 1.1 3.8504 0.9
vary over the length of the reactor, the correlations developed in this study can be applied to these systems only by using average values for the operating parameters. The average operating parameters for an ethylene oxide reactor22 are listed in Table 7. Equations 3-7 were used to calculate the dimensionless operating parameters for the ethylene oxide reactor from the average values. The dimensionless parameters for the ethylene oxide reactor are reported in Table 8. The standard correlation for a nonadiabatic reactor indicates that the dimensionless time to runaway for this system is 0.0145. This dimensionless time corresponds to an actual time to runaway of 87 s. A detailed nonlinear simulation of a fixed-bed ethylene oxide reactor shows that this system has an actual time to runaway of 90 s. The correlation accurately predicts the significant safety concern associated with thermal runaway in ethylene oxide reactors. Although this system has a time to runaway of 90 s, the runaway condition cannot be immediately detected. In industrial applications, ethylene oxide reactor runaway occurs in 15-20 s after detection, making the problem associated with thermal runaway even more difficult. The time to runaway correlations accurately predict the time frame of the time to runaway for the ethylene oxide reactor, indicating that it might be possible to extend the correlations to tubular reactors, although more analysis would be required before making this extension. VIII. Conclusions The work performed in this study focuses on developing an estimate of the time to runaway for an exothermic CSTR. From this estimate of the time to runaway, the severity of the problem associated with thermal runaway can be determined. In this study, correlations for time to runaway in a continuous stirred tank reactor were developed. These correlations relate the time to runaway to the system time constant and system parameters. Separate correlations were developed for adiabatic and nonadiabatic reactors. In the exceptions of operation near the critical heat of reaction parameter or with similar initial reactor and coolant temperatures, correction factors for these correlations were also de-
3730 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
veloped. Judging from several statistical goodness-offit parameters, these correlations predict time to runaway reasonably well for all cases. Judging from the functionality of the correlations developed in this study, several methods of increasing the time to runaway for a given system were identified. These methods of increasing the time to runaway include modifications to the reacting fluid and coolant, as well as to the reactor system itself. Also, the correlations presented in this study were developed assuming a large excess of reactant. Systems without a large excess of reactant will exhibit longer times to runaway. This condition can be achieved by increasing the system pressure or adding an inert component to the reactor feed. Nomenclature A ) reactor surface area for heat transfer (m2) B ) dimensionless heat of reaction BDacritical ) critical dimensionless heat of reaction parameter CA ) reactant concentration (kmol/m3) CA,inlet ) reactant concentration in the inlet stream (kmol/ m3) cp ) heat capacity of the reacting fluid (kJ/kg K) CFNC ) correlation correction factor for near-critical systems CFST ) correlation correction factor for similar-temperature systems Da ) Damko¨hler number k0 ) reaction frequency factor (1/h) k(Tinlet) ) reaction rate constant at the inlet temperature (1/h) n ) reaction order r ) reaction rate (kmol/h‚m3) St ) Stanton number t ) time (h) T ) reactor temperature (K) Tc ) coolant temperature (K) Tinlet ) reactor inlet temperature (K) U ) overall heat-transfer coefficient (kJ/m2‚h‚K) v ) volumetric flow rate of the reacting fluid (m3/h) V ) reactor volume (m3) XA ) reactant conversion γ ) dimensionless activation energy ∆Hrxn ) exothermic heat of reaction (kJ/kmol) Θ ) dimensionless reactor temperature Θ0 ) initial dimensionless reactor temperature Θc ) dimensionless coolant temperature Θinlet ) dimensionless reactor inlet temperature ΘM ) maximum achievable dimensionless temperature for stable systems Θra ) dimensionless temperature at the onset of runaway λ ) system eigenvalue F ) density of reacting fluid (kg/m3) τ ) dimensionless time τA ) dimensionless time to runaway for adiabatic reactors τNA ) dimensionless time to runaway for standard nonadiabatic reactors τNC ) dimensionless time to runaway for near-critical reactors τST ) dimensionless time to runaway for similar-temperature reactors τλ ) system time constant
Literature Cited (1) Semenov, N. N. Zur theorie des verbrennungsprozesses. Z. Phys. 1928, 48, 571-581. (2) Bilous, O.; Amundson, N. R. Chemical Reactor Stability and Sensitivity II. Effect of Parameters on Sensitivity of Empty Tubular Reactors. AIChE J. 1956, 2 (1), 117-126. (3) Morbidelli, M.; Varma, A. Parametric Sensitivity and Runaway in Tubular Reactors. AIChE J. 1982, 28 (5), 705-713. (4) Morbidelli, M.; Varma, A. Parametric Sensitivity and Runaway in Fixed-Bed Catalytic Reactors. Chem. Eng. Sci. 1986, 41 (4), 1063-1071. (5) Morbidelli, M.; Varma, A. A Generalized Criterion for Parametric Sensitivity: Application to Thermal Explosion Theory. Chem. Eng. Sci. 1988, 43 (1), 91-102. (6) van Welsenaere, R. J.; Fromen, G. F. Parametric Sensitivity and Runaway in Fixed Bed Catalytic Reactors. Chem. Eng. Sci. 1970, 25, 1503-1516. (7) McGreavy, C.; Adderley, C. I. Generalized Criteria for Parametric Sensitivity and Temperature Runaway in Catalytic Reactors. Chem. Eng. Sci. 1973, 28, 577-584. (8) Hosten, L. H.; Froment, G. F. Parametric Sensitivity in CoCurrently Cooled Tubular Reactors. Chem. Eng. Sci. 1986, 41 (4), 1073-1080. (9) Henning, G. P.; Perez, G. A. Parametric Sensitivity in FixedBed Catalytic Reactors. Chem. Eng. Sci. 1986, 41 (1), 83-88. (10) Morbidelli, M.; Varma, A. A Generalized Criterion for Parametric Sensitivity: Application to a Pseudohomogeneous Tubular Reactor with Consecutive or Parallel Reactions. Chem. Eng. Sci. 1989, 44 (8), 1675-1696. (11) Soria Lopez, A.; de Lasa, H. I.; Porras, J. A. Parametric Sensitivity of a Fixed-Bed Catalytic Reactor. Chem. Eng. Sci. 1981, 36, 285-291. (12) 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 (5), 1301-1307. (13) Emig, G.; Hofmann, H.; Hoffmann, U.; Fiand, U. Experimental Studies on Runaway of Catalytic Fixed-Bed Reactors (Vinyl Acetate Synthesis). Chem. Eng. Sci. 1980, 35, 249-257. (14) Kimura, S.; Levenspiel, O. Temperature Excursions in Adiabatic Packed Bed Reactors. Ind. End. Chem. Process Des. Dev. 1977, 16 (1), 145-148. (15) Chemburkar, R. M.; Morbidelli, M.; Varma, A. Parametric Sensitivity of a CSTR. Chem. Eng. Sci. 1986, 41 (6), 1647-1654. (16) Hindmarsh, A. C. LSODE and LSODI, Two New Initial Value Ordinary Differential Equation Solvers. ACM-SIGNUM Newsl. 1980, 15 (4), 10-11. (17) Chambre´, P. L. On the Characteristics of a Non-Isothermal Chemical Reactor. Chem. Eng. Sci. 1956, 5, 209-216. (18) Barkelew, C. H. Stability of Chemical Reactors. Chem. Eng. Prog. Symp. Ser. 1959, 55, 37-46. (19) Adler, J.; Enig, J. W. The Critical Conditions in Thermal Explosion Theory with Reactant Consumption. Combust. Flame 1964, 8, 97-103. (20) Thomas, P. H.; Bowes, P. C. Some Aspects of the SelfHeating and Ignition of Cellulosic Materials. Br. J. Appl. Phys. 1961, 12, 222-229. (21) Riggs, J. B. An Introduction to Numerical Methods for Chemical Engineers, 2nd ed.; Texas Tech University Press: Lubbock, TX, 1994. (22) Gudekar, K. G.; Riggs, J. B. Bifurcation and Stability Analysis of an Ethylene Oxide Reactor. Ind. Eng. Chem. Res. 2003, 42 (14), 3285-3293.
Received for review August 11, 2003 Revised manuscript received October 27, 2003 Accepted October 28, 2003 IE0306574