Microsoft Word - 01_R.doc HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 35. pp. 19-30 (2007) RUNAWAY OF CHEMICAL REACTORS: PARAMETRIC SENSITIVITY AND STABILITY F. SZEIFERT, T. CHOVÁN, L. NAGY, J. ABONYI, P. ÁRVA University of Pannonia, Department of Process Engineering, H-8201 Veszprém, P. O. Box 158, HUNGARY Although the runaway phenomenon is well known, an exact consistently used definition does not exist. Present paper is focused on the relation of reactor runaway and parametric sensitivity to stability. The relationship between criteria for reactor runaway and for thermal stability is also pointed out. Based on this critical review new methods based on parametric sensitivity and stability analysis for reactor runaway analysis are proposed in the paper. General runaway criteria based on the application of Ljapunov’s indirect method in the geometric and phase space have been developed. To illustrate the relation of the proposed reactor stability analysis to the commonly used runaway criteria a first order reaction was chosen. The application of the suggested method is also presented on a more complex problem, the Uckron-I test problem. Keywords: runaway, reaction engineering, stability, safety, chemical reactors Introduction Today process control systems allow more and more extensive utilisation of the potential capacity of chemical processes. Optimal operating conditions are pushed to the constraints determined by the laws of physics and chemistry. The most critical safety limit is the runaway condition for reactors with exothermic chemical reactions. Since Semenov’s pioneer work [1] most of the research in reactor analysis is related directly or indirectly to reactor stability and/or reactor runaway. Initially the question was whether reactor runaway is related to reactor stability or it is something completely else. Bilous and Amundson were the first who treat the problem as parametric sensitivity and thermal runaway distinguishing it from classical reactor stability [2]. Schmitz (1975) explained in his review [3] that sensitivity is less precisely defined than stability and the two phenomena might be connected. The first runaway criterion based on the empirical analysis of the temperature profile was suggested by Barkelew [4]. Dente and Collina considered the appearance of an inflection point preceding the temperature maximum (ie. where the second derivative of the temperature with respect to the length of the tubular reactor becomes negative) as a runaway criterion [5]. The same criterion was applied to an industrial problem by Berty et al. four years later [6, 26]. Hlavacek et al. employed an analogous criterion (second time derivative) in the heat explosion theory. Adler and Enig suggested that instead of the length or time domain the analysis be conducted in phase space (temperature versus conversion) and the criterion be the point where the second derivative of the temperature with respect to the conversion becomes negative [7]. Van Welsenaere and Froment called this latter one as “the first criterion” and the previous one as “the second criterion” [8], and found that they are close to each other while the second criterion is a bit more conservative. Applying the method of isoclines Morbidelli and Varma gave the necessary and sufficient conditions [9]. To “measure” the sensitivity Bilous and Amundson introduced the derivative of the maximum temperature along the length (or time) with respect to Semenov’s number [2]. This approach was applied first by Lacey [10] as well as by Boddington et al [11] to determine a runaway criterion. This sensitivity measure gives a maximum curve with respect to the parameter (Semenov’s number). If the parameter is less than the critical value belonging to the maximum then runaway does not occur. Runaway may take place if the parameter is higher. The idea was generalized by Morbidelli and Varma [12]. Introducing the sensitivity Φ * Φ ∂ ∂ T s = and the relative sensitivity Φ * *Φ ∂ ∂ T T S Φ = , where T* stands for the maximum temperature and Φ for an arbitrary input parameter, they defined the runaway criterion based on the maxima of the sΦ – Φ and the SΦ – Φ functions respectively. Their numerical calculations 20 demonstrated that for systems “susceptible to runaway” these maxima practically coincide independently from the selection of Φ and correspond to the value calculated by Adler and Enig’s criterion. At the same time in case of systems less “susceptible to runaway” the maxima can be significantly different or even they can be present when runaway does not occur at all. This observation rises questions about the theoretical soundness of the application of the maxima as runaway criterion in spite of the fact that the authors consider the criteria based on the maxima of sensitivity as important intrinsic attributes of runaway (intrinsic criterion) [12]. Methods for the mainly approximate calculation of various criteria for different reactions and reactors are discussed in several papers, e. g. by Balakotaiah and Kodra [13] as well as Morbidelli and Varma [14]. Asymptotic expressions are often used for the approximation. Methods of catastrophe and singularity theories are more and more often used for studying complex reacting systems. Excellent reviews of these approaches were published by Razón and Schmitz [15], as well as by Doherty and Ottino [16]. Two important sources of the researchers’ motivations are the following: ● foundation of the solution of important industrial problems is expected (reactor design, operation and safety). E.g. in the paper of Luo et al [17] critical ignition, extinction, and transition temperatures and the stable criteria of temperature ● a large variety of processes (with different reactions and different types of reactors) are studied (Zaldivar et al [18]) and they are inherently (very) nonlinear (methods of system analysis can be applied with limitations only). Avoiding runaway is requested from industrial considerations for the sake of “controllability” and for the protection of the catalyst (high amount of heat is released in a short section). The application of runaway criteria in reactor control for building a model-based control strategy is discussed by Szeifert et al [19]. Although the runaway phenomenon is well known, an exact consistently used definition does not exist. Present paper is focused on the relation of reactor runaway and parametric sensitivity to stability. The relationship between criteria for reactor runaway and for thermal stability is also pointed out. In this study Ljapunov’s indirect method that is well known and widely used for dynamic system analysis is applied (Perlmutter, [20]); and it is shown that the commonly used runaway criterion can be justified based on this approach. To allow the comparison of different criteria and to illustrate the application of stability analysis a simple process, a model of a homogenous tubular reactor (or an analogous continuous stirred tank reactor) was chosen. The application of the suggested method is also presented on a more complex problem, known as UCKRON-I test problem (Bashir et al, [21]). Problem definition A commonly used reactor model can be given in the following form: r d cd −= τ (1) ( ) [ ]1,0, ∈−−= ταβ τ w TTr d Td ; (2) where the initial conditions are: τ = 0, c = c0, T = T0 (3) and τ stands for the independent variable which is expressed by the dimensionless length (ℓ · V/F) in case of a steady state tubular reactor or by the time itself in case of a batch reactor. The constitutive equation completing model (1-3) is the following for the rate of reaction: ( ) ( )c T Tcr ϕ δ γ ⋅⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −= exp, . (4) Instead of the reactant concentration c the conversion x(τ) = 1 – c(τ)/c(0) can be used too: ( ) .00 0,= ,0 = = x r d xd c τ τ (5) The explanation of the model parameters is given in the notation list. The runaway phenomenon occurs only in case of exothermic reactions, that is defined by β > 0. Both model (1-4) and model (2-5) describe the concentration (or conversion) and temperature profiles versus the reactor length in case of tubular reactors or versus time in case of batch reactors. A part of the analysis is conducted in the so called phase space. The transformation into the phase space is accomplished by eliminating the independent variable τ (Eq. (2) is divided by Eq. (5)): ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − −≡ −− = r TT r TTr xd Td c ww αβ αβ 0 1 (6) the initial condition is: x = 0, T(x) = T0. The runaway phenomenon is illustrated by Fig. 1. A first order chemical reaction, φ(c) = c was chosen as an example and the following parameters were used: α = 5 l/h, β = 180 m3 K/kmol, γ = 20, δ = 6000 K, c0 = 1 kmol/m 3, T0 = 300 K. (7) (parameters were chosen on the basis of practical experience). Fig. 1 shows that the temperature profile in the steady state tubular reactor (or in a batch reactor) changes significantly for small changes in the operating conditions (eg. changing the value of Tw). 21 260 300 340 380 420 0.0 0.2 0.4 0.6 0.8 1.0 τ [-] T [K ] 270 275 280 285 Tw=290 K 282.4 Figure 1: Parametric sensitivity of a tubular reactor It is well known that model (1-4) gives a unique solution for fixed parameters. That is why Bilous and Amundson distinguished the reactor runaway problem from the multiplicity and classical stability analysis of the continuous tank reactor and interpreted it as parameter sensitivity. For this reason it is practical to determine the parametric sensitivity form of the above model. One among the operating parameters of the model, the Tw external temperature is considered in the further studies. This step does not lessen the generality of the method since the procedure is the same in case of any other parameter. Accordingly the following sensitivities are defined: ( ) , wT c u ∂ ∂ τ = (8) ( ) wT T v ∂ ∂ τ = . (9) Differentiating Eq. (1-2) with respect to Tw the following sensitivity model is obtained: ( )vrur d ud Tc ⋅+⋅−=τ (10) ( ) ( )1−−⋅+⋅= vvrur d vd Tc αβτ ,0 00 = == v uτ (11) where rc and rT are the partial derivatives of the reaction rate with respect to the concentration and the temperature respectively. Eqs (10-11) can be completed with the differential Eqs (1-2) and must be integrated together with those. A differential equation for the calculation of the w(x) = ∂T/∂Tw sensitivity measure in the phase space can be obtained by differentiating Eq. (6) with respect to Tw: ( ) w r Td xd rrTTr rxd wd c xTw 2 0 1 ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⋅+−− −= α α τ = 0, w = 0, (12) where rx stands for the partial derivative of the reaction rate with respect to the conversion ( cx rrc −= 0 1 ). Eq. 12 and Eq. 6 together form a closed system. Stability approach to runaway Let us consider the following system of nonlinear differential equations of the state variables, x(τ): ( )xf d xd = τ . (13) Applying Ljapunov’s indirect method the stability analysis of Eq. (13) is reduced to an eigenvalue analysis of the Jacobian ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = x f J ∂ ∂ for function f(x): │J – λI │ = 0. (14) If the λ1, λ2, … λn roots of Eq. (14) are negative (in case of real roots) or their real components are negative (in case of complex roots) then Eq. (13) is stable at x(τ) [16]. Applying the method on the set of differential Eqs (1-2) the following Jacobian and eigenvalues are obtained respectively: ( ) ( ) ,⎥⎦ ⎤ ⎢⎣ ⎡ − −− = αββ τ Tc Tc rr rr J ( ) ( ) 2 42 2,1 cTcTc rrrrr αβαβαλ −−+±−+− = . (15) 22 Then the conditions for stability can be expressed as: rc + α ≥ β rT (16) If condition (16) is not satisfied then differential Eqs (1-2) are unstable in Ljapunov’s sense. It means also that terms of eλτ, / Re(λ) > 0 appear in the solution of c(τ) and T(τ) instead of eλτ, / Re(λ) < 0 which is the case in the stable region. The Jacobian of the sensitivity model (10-11) is Eq. (15) as well. It means that if Eq. (16) is not satisfied then terms of eλτ, / Re(λ) > 0 are present in the approximate solution for the u and v sensitivities causing significant changes in the profile. Criterion (16) can be “explained” qualitatively too. As exothermic reactions progress, the rate of reaction is generally decreased by the cooling (α) and the “dependence” of the reaction rate on the concentration (these two together form the L. H. S. of Eq. (16)) while it is increased by the “dependence” on the temperature. If the L. H. S of Eq. (16) is the larger then positive feedback is not present; however if the R. H. S is the larger then a positive feedback is formed through the rate of reaction causing instability. Let us apply the stability analysis for model (6) expressed in the phase plane. The stability analysis becomes simpler since the number of state variables is only one. ( ) ( ) .02 < ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ +−− −≡ ≡⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −−≡⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −− r Td xd rrTTr r TT Td d r TTr Td d xTw ww α α αβ (17) After the substitution and rearrangement the following criterion is obtained: ( ) Tww c r TT r TTr rr β β αβ β ≥ − + −− . (18) Obviously this one is different from criterion (16) that is the spatial (length for a tubular reactor, time for a batch reactor) stability criterion and the one expressed in the phase space are not the same. The relation of the two different criteria can be easily determined. For increasing temperature (see balance (2)): βr > α(T – Tw) Considering the above the relationship between the two L. H. S terms of Eq. (16) and Eq. (18): ( ) ccw rr TTr r > −−αβ β α β > − wTT r . (19) Consequently Eq. (16) is always more conservative i. e. the stability in phase space always follows from the spatial stability while inversely does not! Comparing Eq. (12) and Eq. (17) it can be concluded that the relation between sensitivity and stability in the phase space is the same as it is in the spatial analysis. Comparison of runaway criteria The criterion suggested first by Adler and Enig (criterion 1 - van Welsenaere, Froment) can be obtained by differentiating Eq. (6) with respect to x: ( ) 0 1 2 2 0 =⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − −≡⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −− = dx dT r TT dT d r TTr dx d dx Td c ww α αβ (20) Since 0≠ dx dT , Eq. (20) is the same as the result of the stability analysis in the phase space, Eq. (18). It also means that criterion 1, which mainly relies on intuitive bases according to Morbidelli and Varma [12], can be justified theoretically by Ljapunov’s stability analysis. The criterion suggested first by Dente and Collina (criterion 2 - van Welsenaere, Froment) can be obtained by differentiating Eq. (2) with respect to τ: 0 2 2 =−⎟ ⎠ ⎞ ⎜ ⎝ ⎛ += τ α ττ β τ d dT d dT r d dc r d Td Tc . After the substitutions the criterion for avoiding runaway is: ( ) Tw c r TTr rr βα αβ β ≥+ −− . (21) This criterion gives the same result as the “Slope Condition”, dT qd dT qd genrem > (Berty, 1982), which is related to the generated and transferred heat flows, does. Applying the slope condition on the system of Eqs (1-2) the following equivalencies are obtained: ( ) ., ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ −− −→→ w cT genrem TTr r rr dT qd dT qd αβ βα showing that the slope condition differs from the above criterion only in a multiplication with a constant. Berty et al. [22] showed that the criterion, called “Dynamic Condition” by Gilles and Hoffmann [23], can be expressed in the following form, using the original notation: Tc genrem c m T q dT qd ∂ ∂ ∂ ∂ +〉 . Applying the appropriate notation for the system (1-2) the following substitutions are obtained: .,, c T T c genrem r c m r T q dT qd −→→→ ∂ ∂ β ∂ ∂ α Eq. (16) differs from the above criterion only in a multiplication with the same constant value concerning every terms; that is the two criteria is the same! 23 The criterion used at the maximum temperature in the practical design is also included in the study for a more complete comparison: ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ≡≤− T w r rT TT δ 2 . Bashir et al. [21] suggested that this criterion be evaluated at the inflection point instead of the maximum temperature. The above criteria are summarized in Table 1. Criteria 1 and 3 have been compared in the previous section. If the conditions are the same, criterion 2 obviously indicates the runaway always between criteria 1 and 3, i. e. it is more conservative than the first criterion and less conservative than the third one. Criterion 4 used in practical design is especially conservative if the second right hand term is the dominant in criterion 1. Table 1: Runaway criteria for model (1-2) Id Criterion Stability Mathematical form 1 First criterion: inflection in phase space [7] Ljapunov in phase space, Eq. (18) ( )w c w T TTr rr TT r r −− ⋅ + − ≤ αβ 2 Second criterion: inflection in geometric space [7], Eq. (21), “Slope Condition” [22] - ( )w c T TTr rr r −− ⋅ +≤ αββ α 3 “Dynamic Condition” [23] Ljapunov in geometric space (or time), Eq. (16) ββ α c T r r +≤ 4 Practical design (at the hot spot temperature) [21] - w T TT r r − ≤ The relation between the different criteria and the parametric sensitivity measures was also studied. Using the (1-3, 10-11) sensitivity model and assuming that a first order reaction takes place, the u(τ), v(τ) parametric sensitivities can be determined numerically along τ (length or time) and shown on Fig. 2. Naturally the character of the functions changes with the TW parameter. For the sake of obtaining a simpler relationship Morbidelli and Varma evaluated the above sensitivities at the maximum temperature (ST). This sensitivity is shown as a function of TW by the dotted line on Fig. 3. In the literature almost exclusively the hot spot temperature, T* is used for characterizing the reactor operation in relation with thermal runaway. The question is whether the T* value has a distinguished role indeed or some other value, which characterizes the complete reactor regarding runaway and can be measured or calculated easily, can be used too? -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0 0.2 0.4 0.6 0.8 1 τ (length, time) [-] u [k m ol /m 3. K ] -40 -20 0 20 40 60 80 100 120 v [- ] 280285Tw=290 K u v 290 285 280 Tmax (Tw=290 K) . . . . Tmax (Tw=285 K) Figure 2: Sensitivity values versus τ. 24 0 20 40 60 80 100 120 276 278 280 282 284 286 288 Tw [K] S 0 5 10 15 20 25 30 S TLjapunov in geometric space (278.4 K) S ST Inflection in geometric space (281 K) Ljapunov in phase space (282.4 K) . . Figure 3: Parametric sensitivities versus operational parameter Let us express the average rate of reaction in the reactor: ( ) ( )( )∫= 1 0 , τττ dTcrY (22) Based on Eq. (5) the relationship between Y and the conversion can be given as: ( ) Y c x 0 1 1 = . (23) Then a characteristic sensitivity measure can be expressed in the following relative form: w w T Y Y T S ∂ ∂ ⋅= . (24) Applying the (1, 8, 24) relationships S can be calculated from the u(τ) sensitivity: ( ) ( )1 10 u cc T S w ⋅ − −= (25) The solid line on Fig. 3 is the S versus Tw function. Comparing the two sensitivity measures the following conclusions can be drawn: ● the tendency of the two sensitivity measures are similar i. e. the T* value does not have a distinguished role considering runaway; ● the calculation of S value is simpler than that of ST (large numerical error in calculation of T *); ● the maxima of the two functions are close to each other (the higher the inclination for runaway is the closer they are); the difference is theoretically important since it proves that the location of the maximum depends on the choice of the sensitivity measure. To show how the sensitivity measures change for the different criteria summarized in Table 1, the Tw values where the criteria give warning are marked on Fig. 3. The very conservative character of criteria 4 can be noticed here too. In fact the sensitivities start to increase significantly at criterion 3; their values are considerably high at criterion 2; and they reach their maximum values in the narrow neighborhood of criterion 1. The c – T relationships of criteria 1-3 are shown on Fig. 4 for a first order reaction. From Table 1 it is clear that Tw is a parameter for criteria 1 and 2; i. e. each criterion provides a set of curves. On Fig. 4 these curves are given at the critical values of Tw (the operating curves just intersect the criterion curve at the runaway point). Tw is not involved in criterion 3 therefore it gives only a single curve. The operating curve at Tw = 278.4 K does not intersect even the most strict criterion curve 3 so it is considered as runaway less operation. At Tw = 281 K positive exponent components appear in the solution for the state variables along the length (or time), i. e. in Ljapunov’s terms the solution is unstable in space. At the operating point Tw = 282.4 K an inflection is present along the length (or time), however the system is still just stable in phase space. At operating points Tw > 282.4 K the system is in runaway state according to all of the criteria. Fig. 5 shows the same on the T – x plane. From the c – T functions it can be seen that if the c0 initial concentration is less than a critical value, ccr runaway does not occur at all. Fig. 4 and 5 reflect well the relationship between criterion 1-3 too. As it was shown criterion 3 is independent from the Tw value; at the same time as Tw increases, criterion curve 2 approaches curve 3 while curve 1 moves away from it. The very conservative character of criterion 4 has been criticized by several researchers. From criteria 1-3 the third one can be suggested for practical design since 25 it is only slightly more conservative than criteria 1 and 2 and at the same time its preventive character assures reasonable safety. Based on the above example the following conclusions can be drawn: ● Explaining reactor runaway as parametric sensitivity leads to theoretical uncertainty. Although their practical useness is indisputable the criteria obtained this way depend on the selected variable therefore the runaway can not be exactly defined ● The best way is to define the runaway as a stability problem (not in the sense of classical reactor stability). The sudden change in the state variables is then only a consequence of the stability problem. The reason is reaching the stability limits and the result is the change of the shape of functions of state variables with different time delays. (In the case study the stability criterion coincides with the inflection one. This occurs in the case of one independent variable only. In case of more independent variables the criteria do not coincide.) ● The stability analysis must be conducted in the original geometric (or time) space. The indication of runaway may be delayed in transformed spaces. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 270 290 310 330 350 370 T [K] c [k m ol e/ m 3] Tw=277 K 278.4 280 Ljapunov in phase space (282.4) 282.4Ljapunov in geometric space (278.4) 284 Inflexion in geometric space (281) 281 282 Figure 4: Criteria in the c – T plane 260 300 340 380 0.0 0.2 0.4 0.6 0.8 1.0 x (conversion)] [-] T [K ] Tw=277 278.4 280 Ljapunov in phase space (282.4) 282.4 284 Ljapunov in geometric space (278.4) Inflexion in geometric space (281) 281 282 Figure 5: Criteria in the T – x plane 26 Application for a complex system The primary limitation of different runaway criteria published in the literature is that they are generally derived only for simple or considerably simplified (reaction) systems. At the same time the stability analysis described in Section 2 requires only differentiation and calculation of the eigen values. The application of the method therefore does not raise mathematical difficulties. This fact is illustrated using a reactor modeling test problem known as Uckron-I in the literature (Berty [24]). The Uckron-I Test Problem was developed by Berty, Lee, and Szeifert at the Chemical Engineering Department of The University of Akron and by Cropley at the Research and Development of Union Carbide Corporation. Uckron-I is aimed at the computer simulation of a catalytic reaction system for a number of groups to gain experience in chemical reaction engineering. The applied model illustrates some of the complexities of a real system and yet the mathematics involved can be easily handled. The net chemical reaction taking place in the specified industrial methanol synthesis reactor is the following: 2H2 + CO CH3OH . (26) The homogenous one-dimensional axial-mixing-free reactor model for the total mole mass, the components, the enthalpy and the momentum is the following (Berty et al [24]): ( ) rV d cFd 2−= l (27) ( ) { }MCHirV d cFd i i ,,, ∈⋅=ν l (28) ( ) ( )wrp TTAUrVHd Td ccF −−−= Δ l (29) c p f ad F f d pd 2 2 2 ⋅ −= ρ l (30) where ( ) 2.0 3 2.11 8.6 − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⋅ ⋅− = a Fd f p νε ε is the HICKS’ correlation, while ℓ ∈ [0,1]. The rate equation (26) of the net reaction is determined by analysing of the constituent processes and fitting to the data relevant to those processes (Szeifert et al [25]): M C M H pK p p Kp Kr ⋅+ − = 3 2 1 1 , (31) 3,2,1,exp =⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −= i TR E AK iii . The model equations (27-31) can be completed with a suitable state equation (For the sake of simplicity the ideal gas law is applied here). The initial conditions of the system of differential equations are specified at ℓ = 0 in the problem statements (Berty et al [24]): ℓ = 0, F0, ci, 0, T0, p0 are given. (32) Since the number of net reactions is only one the number of state variables F, cH, cC, cM, T, p can be reduced. Since H2 is in excess in the mixture entering the reactor the following form of conversion is introduced: 0,0 0,0 C CC cF cFcF x − = . (33) Based on the (27-28) balance equations and the initial conditions the other state variables can be calculated as functions of the x, p, T state variables: { }MCHic cxc cxc c C Cii i ,,,2 0,0 0,0, ∈ − + = ν (34) . 2 0,0 0 c cxc F F C−= (35) An alternative form of Eq. (34) can be expressed with partial pressures: RT c pp cxc cxc p C Cii i =− + = /, 2 0,0 0,0, ν . (36) Differential Eqs (27-29) can be substituted by the following differential equations applying to x and T: r d dx cC = l 0, (37) ( )[ ], 2 1 wTTrxd dT −− − = αβ κl (38) ℓ = 0, x(0) = 0, T(0) = T0 initial conditions, where ( ) 0, 0 0,0, ,, CpC r pC c c cc H ccV AU = Δ− = ⋅ = κβα . From Eqs (37-38) the differential equation in phase space can be derived: ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − −⋅ − = r TT xdx dT c w C αβ κ 2 11 0, . (39) Model (37-39) differs formally in the factor ( ) ( ) ( )w Cx wC T TTr crr x TT r c r r −− − −+ −⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −≤ αβ κ α 0, 0, / 2 2 1 from model (2-5) of the simple system. This factor expresses the distortion effect due to the change in number of moles in the chemical reaction. According to Ljapunov’s indirect stability criterion the derivative of R. H. S. of Eq. (39) with respect to T must be less than zero. After rearrangement the corresponding runaway criterion is the following: 27 ( ) ( ) ( )w Cx wC T TTr crr x TT r c r r −− − −+ −⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −≤ αβ κ α 0, 0, / 2 2 1 . (40) Comparing the (40) relationship to the first criterion in Table 1 the distortion effect due to the change in number of moles can be described by the following inequalities: 1 2 1 0, <− Cc r α (41) 0 2 1 12 > − <>− κ κ xifx . Consequently one of the two positive terms of R. H. S of Eq. (40) decreases however the other term can increase therefore the relationship between the two terms determines whether the change in number of moles increases or decreases the risk of runaway. The rT and rx derivatives required for the calculation of the criteria can be obtained by differentiating the rate equation of the reaction: 2TR E rrT ⋅= , (42) where the “apparent” energy of activation is: , 1 233 3 1 E p p Bp p p B E pK pK EE C M H C M M M ⋅ − −⋅ + −= ( ) ( ) ⎪⎩ ⎪ ⎨ ⎧ +⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ++− − −= MHX ppKK r pp x r 22 2 1 3 1κ ( ) MC M M C pK K p p p p K x 3 1 2 112 1 +⎪⎭ ⎪ ⎬ ⎫ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + − + (43) According to Eqs (42-43) the (40) criterial equation determines a set of curves with respect to parameter Tw. Any of the curves divides the plane into two regions; one of them is the region where runaway does not occur the other is the runaway region. Fig. 6 shows the criterial curves with respect to parameter Tw in the T – x plane. It can be seen that under a critical value of the conversion the runaway zone for a given conversion is an intermediate temperature region. Both the lower and higher temperatures are outside of the zone. The reason for this is that the reaction is an exothermic reversible one. The criterial curve at Tw = 485.95 K and several operating curves at different values of parameter Tw are given in the cC – T plane on Fig. 7. The criterial curve at Tw = 485.95 K and several operating curves at different values of parameter Tw are given in the cC – T plane on Fig. 7. It can be seen that the operating curve and the criterial curve intersect at Tw = 485.95 K. If Tw is less than this value then the reactor operates without runaway; if Tw is larger then – according to the first criterion – the reactor operates in the runaway region. In case of this system a change of a few tenths K in Tw results in a tremendous change in the reactor temperature (see Fig. 8). The inflection point of the temperature along the length is found at Tw = 485.3 K (criterion 2); positive exponent roots are obtained at Tw = 485.1 K (criterion 3). The specified system is “apt” to runaway that is indicated by the fact that the three criteria give signal in 1 K range. The criterial curves indicating runaways are very important in reactor design. As explained in the discussion of a simple reaction the application of the third criterion is suggested for design since it is only slightly more conservative and with the appearance of the positive exponent roots it practically predicts the imminent changes. 500 520 540 560 580 0.0 0.1 0.2 0.3 0.4 x (conversion) [-] T [K ] runaway zone Tw=475 K 480 485 490 Figure 6: Criterial curves in the T – x plane 28 0.4 0.5 0.6 0.7 0.8 470 490 510 530 550 570 T [K] C c [k m ol /m 3] x=0 Criterion (Tw=485.95 K) Tw=473 K 483 485 485.8 485.9 486 Figure 7: Runaway of the specified reactor in the cC – T plane Reaction (26) is a reversible exothermic one. In reactor design for this type of reactions the determination of the state variables corresponding to the maximum rate of reaction plays an important role too. The rate of reaction for a given composition is between zero (equilibrium) and its maximum value depending on the current value of the temperature. The two extrema can be determined the following way: r = 0 (equilibrium) (44) rT = 0 (max. reaction rate - necessary condition). (45) Eqs (44) and (45) give the corresponding relationships in the T – x and ci – T planes according to Eqs (41) and (42) respectively. Each curve can be constructed by solving the corresponding nonlinear equation. Fig. 9 shows the runaway curve corresponding to the third criterion as well as the equilibrium and maximum curves of the rate of reaction in cC – T diagram. It designates the main limits for the design. In case of low conversion the temperature providing higher rate of reaction is restricted by the risk of runaway. Over a critical xcr conversion runaway does not occur therefore the temperature can be increased to approach the one corresponding to the maximum rate of reaction. The approximation of the “ideal operating curve” shown on Fig. 9 is a part of the design process and it is outside of the scope of present study. 460 500 540 580 620 0 2 4 6 8 10 12 14 Length [m] T [K ] Tw=503 K 486 485.9 483 473 Figure 8: Temperature gradients along the reactor length 29 0.3 0.4 0.5 0.6 0.7 0.8 490 510 530 550 570 T [K] C C [k m ol /m 3] runaway r=0 rmax x=xcr x=0 x=0 Figure 9: Design diagram Conclusions The most important criteria for reactor runaway are summarized in the paper. It is proved that reactor runaway approached as parametric sensitivity can be interpreted as a consequence of stability problems. After the reactor gets into the unstable region of the geometric (length) space (criterion 3) – as a consequence – a convex temperature profile appears soon (criterion 2). “Going on” into the unstable region the phase space instability criterion (criterion 1) is reached shortly. The absolute values of the different parametric sensitivities reach their maximum values in this narrow range; however from the fact itself that the parametric sensitivity is at its maximum reactor runaway does not follow. The correspondence and relations of reactor runaway, stability and thermal stability criteria are illustrated on the example of a simple reactor. It is proved that the maximum temperature (hot spot) does not have a distinguished role in the analysis of runaway. This also comes from the fact that the component and enthalpy balances form a mutually linked system. The application of Ljapunov’s stability criteria for reactor runaway is demonstrated also on a more complex problem showing that complexity does not restrict the use of the method. To avoid runaway the application of the third criterion is suggested for reactor design or control since it is only slightly conservative (compared to the first and second criteria) and it is the first among criteria 1-3 to predict the appearance of dangerous operating conditions. NOTATION A heat transfer area, m2 a cross section, m2 c concentration, kmol/m3 cp specific heat, kJ/m 3K dp particle diameter, m E energy of activation, kJ/kmol F volumetric flow rate, m3/h fc friction factor (–ΔHr) heat of reaction, kJ/kmol J Jacobian matrix K kinetic parameters ℓ dimensionless length, (-) m mass balance function, kmol/m3h p pressure, bar genq generated heat flow, kJ/h remq transferred heat flow, kJ/h R gas constant, kJ/kmol K r rate of reaction, kmol/m3h rc, rT, rx derivative of reaction rate with respect to c, T and x respectively S, ST relative sensitivity s absolute sensitivity T temperature, K T* hot spot temperature, K Tw wall temperature, K u,v parametric sensitivities x conversion Y average rate of reaction, kmol/m3h ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = pcV UA ρ α parameter, 1/h ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − = p r c H ρ β Δ parameter, m3K/kmol γ kinetic parameter, (-) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = R E δ kinetic parameter, K ε void fraction, m3/m3 ϕ(c) concentration function of the rate of reaction, kmol/m3h ν dynamic viscosity, kg/m s 30 νi stochiometric coefficient ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = 0, 0 cc c κ initial ratio of concentrations, (-) ρ density, kg/m3 SUBSCRIPTS C CO H H2 M CH3OH 0 inlet cr critical REFERENCES 1. SEMENOV N. N.: Zur Theorie des Verbreennungs- prozesses, 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, AIChEJ., 1956, 2, 117-126 3. SCHMITZ R. A.: Multiplicity, Stability and Sensitivity of States in Chemically Reacting Systems, A Review. Adv. Chem. Ser., 1975, 148-156 4. BARKELEW C. H.: Stability of Chemical Reactors, Chem. Eng. Prog. Symp. Ser., 1959, 25, 37-46. 5. DENTE M., COLLINA A.: Il Comportamento dei Reattori Chimici a Flusso Longitudinale nei Rigvardi della Sensitívitá, Chim. e Industria, 1964, 46, 752- 761 6. BERTY J. M., BRICKER J. H., HAMBRICK J. O.: Parametric Sensitivity and Stability of Staged Adiabatic Reactors with Interstage Coolers, Symposion on Stability and Control of Reaction Systems: Part III, Preprint 10E, St. Louis, Missouri, Feb. 1968 7. ADLER J., ENIG J. W.: The Critical Conditions in Thermal Explosion Theory with Reactant Consumption, Comb. Flame, 1964, 8, 97-103 8. VAN WELSENAERE R. J., FROMENT G. F.: Parametric Sensitivity and Runaway in Fixed Bed Catalytic Reactors, Chem. Eng. Sci., 1970, 25, 1503-1516 9. MORBIDELLI M., VARMA A.: Parametric Sensitivity and Runaway in Tubular Reactors, AIChEJ., 1982, 28, 705-712 10. LACEY A. A.: Critical Behavior for Homogeneous Reacting Systems with Large Activation Energy, Int. J. Engng. Sci., 1983, 21, 501-515 11. BODDINGTON T., GRAY P., KORDILEWSKI W., SCOTT S. K.: Thermal Explosions with Extensive Reactant Consumption: A New Criterion for Criticality, Proc. R. Soc., 1983, A 390, 13-30 12. MORBIDELLI M., VARMA A.: A Generalized Criterion for Parametric Sensitivity: Application to Thermal Explosion Theory, Chem. Eng. Sci., 1988, 43, 91-102 13. BALAKOTAIAH V., KODRA D.: Stability Criteria for Chemical Reactors, DYCORD+`92, Preprints, Maryland, USA, 1992, 83-92 14. MORBIDELLI M., VARMA A.: On Parametric Sensitivity and Runaway Criteria of Pseudohomogeneous Tubular Reactors, Chem. Eng. Sci., 1985, 40, 2165-2168 15. RAZÓN L. F., SCHMITZ R. A.: Multiplicities and Instabilities in Chemically Reacting Systems - A Review, Chem. Eng. Sci., 1987, 42, 1005-1047 16. DOHERTY M. F., OTTINO J. M.: Chaos in Deterministic Systems: Strange Attractors, Turbulence and Applications in Chemical Engineering, Chem. Eng. Sci., 1988, 43, 139-183 17. LUO K. M., LU K. T., HU K. H.: The critical condition and stability of exothermic chemical reaction in a non-isothermal reactor, J. Loss Prevention in the Proc. Ind., 1997, 10, 3, 141-150 18. ZALDIVAR J. M., CANO J., ALOS M. A.: A general criterion to define runaway limits in chemical reactors. Journal of Loss Prevention, 2003, 16, 187- 200 19. SZEIFERT F., CHOVAN T., NAGY L.: Process dynamics and temperature control of fed-batch reactors, Computers & Chemical Engineering, 1995, 19(11), 447-452 20. PERLMUTTER D. D.: Stability of Chemical Reactors, Prentice-Hall, Englewood Cliffs, 1972 21. BASHIR S., CHOVAN T., MASRI B. J., MUKHERJEE A., PANT A., SEN S., VIJAYARAGHAVAN, BERTY J. M.: Thermal Runaway Limit of Tubular Reactors, Defined at the Inflection Point of the Temperature Profile, Ind. Eng. Chem. Res., 1992, 31, 2164-2171 22. BERTY J. M., LENCZYK J. P., SHAH S. M.: (1982) Experimental measurement of the Thermal Stability Criteria for the Low Pressure Methanol Synthesis, AIChE Journal, 28, 914-922. 23. GILLES E. D., HOFFMANN H.: An Analysis of Chemical Reactor Stability and Control, Chem. Eng. Sci., 1961, 5, 328 24. BERTY J. M., LEE S. G., SZEIFERT F., CROPLEY J. B.: The “UCKRON-1” Test Problem for Reaction Engineering Modeling, Chem. Eng. Comm., 1989, 76, 9-33 25. SZEIFERT F., ÁRVA P., NAGY D.: Effects of Pore Diffusion on the Synthesis of Methanol, Chem. Eng. Comm., 1989, 76, 157-167 26. BERTY J. M.: Experiments in Catalytic Reaction Engineering, Elsevier, Amsterdam, 1999