- Email: [email protected]
Study of gas–liquid stirred reactor stability taking into account the temperature influence on gas solubility
Chemical Engineering Science 59 (2004) 4137 – 4147 www.elsevier.com/locate/ces
Study of gas–liquid stirred reactor stability taking into account the temperature influence on gas solubility E.F. Stefoglo∗ , I.V. Kuchin, A.V. Kravtsov Institute of Coal and Coal Chemistry, Rukavishnikov, 21, 650610, Kemerovo, Russia Received 3 December 2003; received in revised form 6 April 2004; accepted 25 April 2004
Abstract The paper is devoted to the dynamic behavior and stability of gas–liquid stirred reactor taking into account the temperature influence on gas solubility. Since the rate of gas–liquid processes is very sensitive to concentration of gas reactant dissolved in liquid, even weak fluctuations of temperature can significantly influence on process pass. There are two cases of temperature influence on gas solubility are possible: (1) the solubility decreases with increasing temperature; (2) the solubility grows with increasing temperature. The first case is typical for majority of gases. The second case occurs more rarely but has a great practical importance. It takes place, for example, for the hydrogenation of many compounds in organic solvents (such as benzene, toluene, isopropyl alcohol and others). A model of gas–liquid process has been developed to demonstrate the stability of gas–liquid reactor. It has been shown that the gas solubility behavior has an influence on the form of heat production curve and therefore on the multiplicity of the steady states. The areas of multiplicity and limit cycles were found and the phenomenon of hysteresis in the reactor was shown. A criterion to determine whether the multiplicity is possible under the given conditions was found. By means of an analysis of a mathematical model the stability of steady states of the reactor was studied. 䉷 2004 Elsevier Ltd. All rights reserved. Keywords: Multiphase reactors; Dynamic simulation; Stability; Mathematical modeling; Multiplicity; Absorption
1. Introduction Multiplicity, stability and dynamic behavior of singlephase processes have been well studied in the literature. During recent decades, many papers have been devoted to this subject. The main results in this area are shown in papers (Bilous and Amundson, 1955; Aris and Amundson, 1958; Slinko and Muler, 1961; Frank-Kamenetski, 1969; Olsen and Epstein, 1993 and others) in detail. Prediction of the dynamic behavior of multiphase reactors is usually more complex and less investigated because these systems involve: (1) more than two component balances, i.e. ODEs, (2) mass transfer between the gas and
∗ Corresponding author. Tel.: 7-3842-36-55-61; fax: 7-3842-21-18-38.
E-mail address: [email protected] (E.F. Stefoglo) 0009-2509/$ - see front matter 䉷 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2004.04.041
the liquid phase and (3) gas reactant that comes into the reaction zone independently of the liquid reactant. The existing investigations (Hoffman et al., 1975; Sharma et al., 1976; Beskov et al., 1979; Huang and Varma, 1981; Singh and Shah, 1982; Vleeschhouwer et al., 1992; van Elk et al., 1999) do not take into account that the changing temperature leads to the changes in the gas solubility in the liquid phase. However, very often the real processes are carried out at the unsteady temperature. Changes in the temperature lead to the changes in gas solubility in liquid phase. Dependence of gas solubility on the temperature really can be very significant. The temperature change in 10 ◦ only can lead sometimes to the change of gas solubility several times (such systems as hydrogen–acetone, hydrogen–benzene, etc.). Since the concentration of gas reactant dissolved in liquid is, as a rule, much lower than the liquid reactant concentration, the gas concentration can influence the reaction rate to a great
4138
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
extent. Thus, in the study of thermal stability of a reactor it is very important to take into account the effect of temperature on gas reactant content in the liquid. The dynamics of gas–liquid process are considered in this paper by means of a model accounting for the temperature influence on gas solubility. First-order reaction (r = kC GL ) is considered of the type G(gas) + L(liquid reactant) → P (product) carried out in continuous stirred tank reactor (CSTR).
2. Mathematical model of CSTR taking into account temperature influence on gas solubility For steady-state operating gas–liquid reactor the material and heat balance are: for a gas reactant: ∗ VL (CG − CGL ) − VL 1 kC GL − qL CGL = 0,
(1)
Two cases of temperature influence on gas solubility are possible: 1. solubility decreases as temperature increases. 2. solubility increases as temperature increases. The first case is typical for majority of gases. The Henry’s law constant depends on temperature according to equation H e = H e∞ exp(−Q/RT ). The second case is not so common but is of great practical importance. It takes place, for example, with hydrogen dissolution in organic solvents (benzene, toluene, isopropyl alcohol and others). In this case, the Henry’s law constant must decrease with increasing temperature. This is possible if the heat of absorption Q in expression H e = H e∞ exp(−Q/RT ) is formally considered as a negative magnitude. In the quite wide range of temperature far from the boiling point of solvent the absorption heat Q may be accepted as a constant magnitude not depending on temperature. Introducing the Stanton number as
for a liquid reactant: qL CL0 − qL CL − VL 2 kC GL = 0,
(2)
(3)
where Qin L = qL L T0 CT L is the heat entered the reactor by the liquid flow, [W], Qout L = qL L CT L T the heat removal from the reactor with a liquid flow, [W], QR = VL H r is the heat production as a result of chemical reaction, [W], Kc Fc (T − Tc ) the heat removal from the reactor through heat exchange surface, [W]. ∗ = P /H e), From Eqs. (1)–(2) and Henry’s law (CG Eq. (3) takes the form: Kc Fc Tad · = (T − T0 ) + (T − Tc ). L CT L qL
Kc Fc , L C T L q L
we write the right-hand side of Eq. (4) as
heat balance equation: out Qin L − QL + QR − Kc Fc (T − Tc ) = 0,
St =
(4)
(1 + St)(T − TR ), where TR =
T0 + St · Tc . 1 + St
(5)
3. Results of model solution. Multiplicity of steady states Heat production curve Q1 and heat removal line Q2 are shown in Fig. 1 for the case when gas solubility is
The left-hand side of Eq. (4) represents the heat production (Q1 ) as a result of chemical reaction, while the righthand side represents the heat removal (Q2 ) with outgoing liquid flow and heat exchange. In Eq. (4) the following magnitudes are used:
Tad =
H C L0 — adiabatic heating, L C T L
k∞ e−E/RT · P 2 CL0 H e∞ e−Q/RT ( + k∞ e−E/RT + 1) — conversion of liquid reactant,
=
is the residence time of the solution in the reactor, [s], H e =H e∞ exp(−Q/RT ) is the dependence of Henry’s law constant on the temperature; Q is the heat of absorption, [J/kmol].
Fig. 1. Heat production curves and heat removal line, steady-state multiplicity in CSTR for a the first-order reaction. The gas solubility is reduced as temperature increases.
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
Fig. 2. Heat production curve and heat removal line for the first-order reaction in CSTR. The gas solubility increases as temperature increases.
reduced with temperature increasing. The plot is drawn at the following values of process parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 1 × 106 s−1 ; H e∞ = 1 × 107 Pa m3 /kmol; Q=3×106 J/kmol; E=4.9×107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); = 0.07 s−1 ; P = 101000 Pa; qL = 8 × 10−6 m3 /s; = 1.25 × 103 s; L = 1300 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.
Fig. 3. Areas of steady states multiplicity in gas–liquid reactor in terms of temperatures T and TR = (Q 1 · T − Q1 )/Q . First-order reaction. The gas solubility is reduced as temperature increases. Calculations were made on data of Fig. 1.
At the inflection point (TF in Figs. 1 and 2) a slope of tangent to heat production curve is maximum. If the heat removal line has a slope greater than the maximum slope, then only one intersection point of heat production curve with heat removal line is possible. The multiplicity is not possible in such conditions. The tangency point TF with the maximum slope of tangent is a trifurcation point. Using expression (4) we find the slope of tangent to heat production curve at trifurcation point:
E−Q · P 2 k∞ e−E/RT max RT ( + k∞ e−E/RT max + 1) − k∞ e−E/RT max RTE2 2 1 + St d(Q1 )max max max = , = Tad dT H e∞ e−Q/RT max CL0 ( + k∞ e−E/RT max + 1)2 It has been seen from Fig. 1, that at T > 430 K the heat production is reduced because the gas solubility decreases with temperature increasing. Two positions of heat removal line are shown in Fig. 1a and b—the tangents to heat production curve. Both lines have the same slope and intersect an abscissa axis at points TRa and TRb . Three areas can be indicated in Fig. 1: (1) TR < TRa : in this area only one stable steady state exists with low value of liquid reactant conversion. (2) TRa < TR < TRb : in this area three steady-states occur—two stable states with low and high conversion, one unstable state with middle conversion. (3) TR > TRb : in this area only one stable steady state with high conversion exists. Heat production curve Q1 and heat removal line Q2 are shown in Fig. 2 for the case when with temperature increasing the gas solubility increases (Q < 0). The plot is drawn at the following values of parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 2.2 × 106 s−1 ; H e∞ = 3.5 × 105 Pa m3 /kmol; Q = −2.5 × 106 J/kmol; E = 5 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); = 0.2 s−1 ; P = 101000 Pa; qL = 2 × 10−4 m3 /s; = 50 s; L = 1300 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 275 K; KC = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.
4139
(6)
where Tmax is the temperature corresponding to trifurcation point TF. To find the temperature Tmax corresponding to trifurcation point, it is necessary to solve numerically the following equation: d2 (Q1 ) = 0. (7) dT 2 For data of Fig. 1 Tmax is equal to 348.4 K. Thus, a criterion of a single stable steady-state existence is 1 + St d(Q1 )max > , Tad dT
(8)
1 )max is calculated using Eq. (6). where d(QdT It should be noted, that value 1/Q 1 = 1/(dQ1 /dT ) is the maximum increase in reactor temperature if full conversion of the liquid reactant takes place under the non-adiabatic conditions. The areas where multiplicity can take place are shown in Fig. 3 in coordinates (TR ; 1/Q 1 ) and (T; 1/Q 1 ). As can be seen from Fig. 3, a single, stable steadystate takes place at 1/Q 1 < 0.3 that corresponds to values TR = 311 K, T = Tmax = 348.4 K. From Fig. 3, it follows that at a value 1/Q 1 = 0.7 the multiplicity takes place within the temperature intervals 227.8 < TR < 291.2 K and 310 < T < 387.7 K.
4140
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
Fig. 5. The hysteresis in CSTR for the case when gas solubility decreases as temperature increases (Q > 0). Calculations were made on data of Fig. 1. Fig. 4. Areas of multiplicity in gas–liquid reactor in terms of temperatures T and TR = (Q 1 · T − Q1 )/Q . First-order reaction. The gas solubility increases as temperature increases. Calculations were made on data of Fig 2.
Fig. 4 is drawn similarly to Fig. 3 but for the case in which the gas solubility increases with increasing temperature. We have one stable steady state in this case at 1/Q 1 < 0.5 that corresponds to values TR = 331.1 K, T = Tmax = 371.7 K.
4. Influence of gas solubility on the hysteresis in CSTR Let us consider the reactor operating in a stable steadystate with low value of liquid reactant conversion (Fig. 1). As the temperature T0 of the liquid flowing to the reactor inlet, the heat removal line moves in parallel to itself to the right and at the certain moment will coincide with a tangent to heat production curve (line b in Fig. 1). As T0 increases further, only one steady-state characterized by the high conversion will take place. Thus, at a given value of T0 called the ignition temperature (T0 )ign , a strong increase in the reaction rate has occurred and the steady state has moved to the upper position. For the opposite case, i.e., for a reduction in the temperature T0 of liquid flowing to the reactor inlet, the process occurs at a high conversion until the reaction is extinguished at a certain temperature (T0 )ext . The hysteresis phenomena means that the ignition and extinction temperatures are dissimilar: (T0 )ign = (T0 )ext . The hysteresis diagram (Fig. 5) is calculated from the data of Fig. 1. The hysteresis interval in this case is equal: (T0 )ign − (T0 )ext = 67.9 ◦ C. A similar diagram is obtained for the case in which the gas solubility increases with temperature (Q < 0). The results of this calculation are presented in Fig. 6. Calculations were made on the data of Fig. 2. In this case the hysteresis interval is constant: (T0 )ign − (T0 )ext = 49.5 ◦ C.
Fig. 6. The hysteresis in CSTR for case when with gas solubility increases as temperature increases (Q > 0). Calculations were made on data of Fig. 2.
5. Criterion of steady-states multiplicity in the reactor We now develop a criterion indicating the possibility of multiplicity in gas–liquid reactor, using the approach proposed by Slinko (1961). From Eq. (4) we have heat production (HPR): Tad · , heat removal (HWR): (T − T0 ) +
Kc F c (T − Tc ). L CT L qL
In the steady-state HWR = HPR or HWR/HPR = 1. Let us introduce the following dimensionless values:
= E/(RT 20 )(T − T0 );
T =
RT 20 + T0 ; E
RT 0 = ; E L ; exp −
L = Q/E; S = P /(H e0 CL0 ); H e0 = H e∞ k = k0 exp ; k0 = k∞ exp(−1/). 1 +
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
4141
The differentiation of function F () gives dF d A B + (−q)2 k0 exp 1+ − B(1 − L) (1+) , = (1−L) 2 S k0 exp 1+
F () ≡
(11) where A = (1 + U k 0 · ), B = + k0 exp Fig. 7. The function F () at various values of parameter k0 . Calculations were made on data of Fig. 1.
Then Eq. (4) can be rewritten in the form: E HWR Tad 2 HPR RT 0 Kc Fc E ) + (T − T ) (T − T 0 c 2 C q RT L TL L . = 0
≡F (),
(13)
k0 exp
− B(1 − L) 1 + (14)
By substituting Eq. (13) in Eq. (14), after some transformations we discover: k0 exp 1 + ( − q) =− + k exp + 1 0 1 + (1 + )2
×L − − 1 − − 1.
(15)
(10) U ∗ k0 c 1+U ∗ k0 ;
c = E/(RT 20 ) · (Tc − T0 ); U ∗ = Kc Fc 2 k0 L CT L qL ; ad = Tad E/(RT 0 ). At steady states, the relation HWR/HPR = 1. The intersection points of curve F () and line ad correspond to the steady states of the reactor. We examine the function F () to define when the multiplicity is possible. In Fig. 7 a plot of the function F () at various values of parameter k0 is shown. For different values of k0 there can exist three solutions or just one, depending on the value ad . In Fig. 7 ad = 20.7, and only at k0 = 2.6 three intersection points of function F () with ad are possible. At other values of parameter k0 only one intersection point exists. With an increase in the parameter k0 the extrema of function F () become less expressed and disappear completely at certain k0 value. For such a form of the curve F () only one intersection point of function F () and ad can occur. This means that there is a single steady state for any values of ad , and thus multiplicity is impossible. Therefore, it is necessary to formulate a condition when the function F () loses the extrema. where q =
( − q)
(1 + )2 =0.
If the left-hand side of Eq. (9) note as ad and the above dimensionless complexes introduce into the right-hand side, Eq. (9) takes the form: ( − q)(1 + U ∗ k0 ) · + k0 exp 1+ + 1 ad = (1−L) k0 exp 1+ 2 · S
+ 1. 1 +
To find the extrema of function F () we equate F () to zero. Taking into account that A = 0 and the denominator of expression (11) is not zero, we find: B+
(9)
(12)
The left-hand side of expression (15) is some function of (we note it as f1 ()). The right-hand side of this expression we note as f2 (). Then expression (15) takes the following form: f1 () = f2 (). If f1 () and f2 () intersect each other at some points 1 and 2 , then derivative F is equal to zero at these points, therefore the function F has extrema. If f1 () = f2 (), the functions f1 and f2 do not intersect each other, F = 0 and function F has no extrema, therefore the multiplicity is impossible. Thus, the problem of determination of function F () extrema reduces to the finding of intersection points of the curves f1 () and f2 (). To define whether the intersection of functions f1 and f2 takes place, the following method can be used. We draw two parallel tangents to functions f1 and f2 by such way that the line passing through tangency points 1 and 2 to be perpendicular to the tangents (Fig. 8). Values 1 and 2
4142
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
characterized by the high value of conversion. We note this steady state by point S. Using the expansion proposed by Frank–Kamenetski (1969) we find:
Fig. 8. Functions f1 and f2 dependences on .
corresponding to points 1 and 2 can be found by numerical solution of the following equations set: f1 (1 ) = f2 (2 ), f2 (2 ) − f1 (1 ) = −1. f1 (1 ) 2 − 1
(16)
The first equation means that the tangents to functions f1 and f2 at points 1 and 2 are parallel (their derivatives are equal). The second equation means that the line passing through tangency points 1 and 2 is perpendicular to tangents (if lines are perpendicular, the product of their angular factors is equal to −1). Calculation of 1 and 2 allows to we find values of functions f1 (1 ) and f2 (2 ). The condition when the functions f1 and f2 are not intersected, and therefore the multiplicity is impossible is f1 (1 ) > f2 (2 )
or
2 > 1 .
(17)
6. Stability of steady states of reactor In practice, we deal with non-stationary operation of continuous reactor because the various deviations of process parameters are possible. Therefore it is necessary to know when these deviations do not lead to the change of process mode, i.e. to know the conditions of stable reactor operation. Under the non-stationary conditions the mass balance equation for gas reactant is: dCL ∗ − CGL ) − 1 k CGL − CGL , = (CG dt
(18)
where t = ϑ/ is the dimensionless time; ϑ the current time, ∗ = P /H e the equilibrium [s]; the residence time, [s]; CG concentration of gas dissolved in liquid. On the basis of expression (4) we write an equation for temperature in reactor: dT = (T0 − T ) − St(T − Tc ) + Tad · . dt
(19)
As described above (Figs. 1 and 2), there are no more than three intersection points of heat production curve and heat removal line. One point corresponds to an unstable steady state and two points—to stable ones. The upper stable steady-state represents the special interest because it is
k = kS exp(),
(20)
H e = H eS exp( · L),
(21)
where H es is the Henry’s law constant at the temperature, Ts of the upper stable steady state, [Pa m3 /kmol]; the dimensionless deviation from temperature Ts , = E/(RT 2s )(T − Ts ); L = Q/E the dimensionless parameter; Q the heat of absorption, [J/kmol]; E the activation energy, [J/kmol]; ks the reaction rate constant at the temperature Ts , [s−1 ]; By dividing of Eq. (18) on Eq. (19), using dependences (20), (21) and introducing dimensionless parameters, we find H e Pe·L −CGL −CGL (v1 ks e +1) dCGL s = , d −(1 + St)+ad v2 ksve1 CLCGL −ad s
(22)
0
where ad = Tad E/(RT 2S ) is the dimensionless parameter, s the conversion of liquid reactant at the temperature Ts ; St = Kc Fc /(L CT L qL ) the Stanton’s criterion; 1 , 2 the stoichiometric coefficients. From Eq. (22) the dependence of gas concentration CGL on the dimensionless temperature deviation from a steady state was found. Knowing the dependence CGL on , it is easy to find the function on with the aid of relation:
= ks e · CGL /CL0 .
(23)
In Fig. 9 the dependence on in combination with heat production and heat removal lines are shown. If as a result of some disturbance, the process deviates from stable state (point S) to some state (point A), then the process will return to S by the trajectories presented in Fig. 9. A similar plot for the case when the solubility of gas increases with temperature increasing (Q < 0) is shown in Fig. 10. Solution of Eq. (22) is the integral curves of the following system: dCGL P = − CGL − CGL (v1 ks e + 1) dt H eS e·L ≡1 (CGL , ), v2 ks e CGL d =−(1 + St) + ad − ad s dt v1 CL0 ≡2 (CGL , ).
(24)
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
4143
Fig. 11. Areas corresponding to various types of equilibrium position.
= 0 in the upper steady state, the coefficients of the linearized equations willbe equal: Fig. 9. The dependence and CGL on for the first-order reaction in gas–liquid CSTR. The solubility of gas increases as temperature decreases (Q > 0). Calculations were made on data of Fig. 1.
a11 =
*1 = − − v1 ks − 1, *CGL
a12 =
*1 P =− L − v1 ks · CGLs , H es *
a21 =
*2 ad v2 ks = , v1 CL0 *CGL
a22 =
*2 ad v2 ks · CGLs = −(1 + St) + . v1 CL0 *
The type of equilibrium position is determined by values of known criterions: = −(a11 + a22 ), = a11 a22 − a12 a21 . (26)
Fig. 10. The dependence and CGL on for the first-order reaction in gas–liquid CSTR. The solubility of gas increases as temperature increases (Q < 0). Calculations were made on data of Fig. 2.
The points of phase plane corresponding to equilibrium position (steady state of the reactor) satisfy to the conditions: dCGL d = = 0. dt dt Then
1 (CGLS , S ) = 2 (CGLS , S ) = 0,
(25)
where CGLS , S are values of gas concentration and temperature deviation under equilibrium condition (at the temperature TS ). To study the process stability it is convenient to linearize Eqs. (24), near the temperature TS . Taking into account that
The positive values of parameters and correspond to stable equilibrium positions (Fig. 11). Then the conditions of stability in terms of expressions for coefficients a11 , a12 , a21 , a22 , are: ad v2 ks · CGLs (− − v1 ks − 1) · −(1 + St) + v1 CL0 P ad v2 ks − L − v1 ks · CGLs · > 0, (27) H es v1 CL0 ad v2 ks · CGLs ( + v1 ks + 1) − −(1 + St) + v1 CL0 >0 (28) or using expression (23) and taking into account that S =0:
E v 2 S (− − v1 ks − 1) · −(1 + St) + RT 2S v1 P E Tad v2 kS − L − v1 S CL0 · > 0, (27a) H es RT 2S v1 CL0
4144
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
Fig. 12. Limit cycle, dependence on for the first-order reaction in gas–liquid CSTR.
E v2 S ( + v1 ks + 1) − −(1 + St) + RT 2S v1
> 0. (28a)
The areas corresponding to various types of equilibrium positions in coordinates , are shown in Fig. 11. The point A corresponds to values and calculated on data of Fig. 1, the point B—on data of Fig. 2. Eq. (27) represents a static stability condition, and Eq. (28), the dynamic stability condition. If conditions (27) and (28) are met, the process is stable in steady state. If condition (27) is met and condition (28) is not met, the reactor is dynamically unstable and operates cyclically in the limit cycle around the point S corresponding to the upper steady state (Fig. 12). The plot is drawn at following values of parameters: VL =0.01 m3 ; H =1.7×108 J/kmol; k∞ =1×106 s−1 ; H e∞ = 2.5 × 107 Pa m3 /kmol; Q = 2 × 107 J/kmol; E = 4.9 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); =0.05 s−1 ; P =101000 Pa; qL =9×10−4 m3 /s; =11.1 s; L = 1200 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.The transition to a limit cycle will take place as a result of process deviation from a steady state to point A. The change of process variables (gas dissolved concentration, conversion of liquid reactant, temperature) in the reactor has a cyclical character with a certain frequency and amplitude. The oscillations of values CGL and corresponding to the case of Fig. 12 are shown as a function of time in Fig. 13. We draw the diagram of stability of the given reactor. From expression (4) we find a relation between S and TS :
S =
1 + St (TS − TR ), Tad
where S =
k∞ e−E/RT S ·P 2 −Q/RT S (+k∞ e−E/RT S +1) CL0 H e∞ e
(29)
Fig. 14. Stability diagram for CSTR. Gas solubility increases as temperature increases (Q < 0). (1) zero conversion; (2) full conversion; (3) heat production; (4) area of multiplicity; (5) area of limit cycles.
Dependence (30) is shown as a curve 3 in Fig. 14. Calculations were made at following values of parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 9 × 106 s−1 ; H e∞ = 5 × 105 Pa m3 /kmol; Q = 9 × 106 J/kmol; E = 4.9 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); =0.05 s−1 ; P =101000 Pa; qL =9×10−4 m3 /s; =11.1 s; L = 1200 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K. At zero conversion (S = 0) we have from Eq. (30): TS − TR = 0.
(30)
(31)
This condition corresponds to line 1 in Fig. 14. Line 2 corresponds to full conversion S = 1: TS − TR =
.
Rewriting Eq. (29), we find:
Tad TS − TR = . 1 + St S
Fig. 13. Oscillations of values CGL and during time.
Tad . 1 + St
(32)
Eq. (27a) is a condition of static stability. To find the boundary of static stability region, we equate the lefthand side of Eq. (27a) to zero, substituting instead of S
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
4145
CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.
7. Conclusions Since dissolved gas concentration is a magnitude of several orders lower than liquid reactant concentration, the rate of gas–liquid processes is very sensitive to concentration of gas reactant dissolved in liquid, and so even weak fluctuations of temperature can significantly influence on process pass. In the given paper the dynamics of gas–liquid process in the CSTR with taking into account the temperature influence on gas solubility has been considered for two cases: Fig. 15. Stability diagram for CSTR. Gas solubility reduces as temperature increases (Q > 0). (1) zero conversion; (2) full conversion; (3) heat production; (4) area of multiplicity; (5) area of limit cycles.
the right-hand side of Eq. (29): (− − v1 ks − 1) E v2 (1 + St)(Ts − TR ) × −(1 + St) + v1 RT 2s P 1 + St L − CL0 (Ts − TR ) − − H es Tad E Tad v2 ks × = 0. RT 2s v1 CL0
(33)
Using Eq. (33) we draw a curve 4 in Fig.14, which covers the statically unstable area of reactor operation (area of multiplicity). As a result of similar operations with condition (28a) of dynamic stability, we obtain E ( + v1 ks + 1) − −(1 + St) + RT 2s v2 (1 + St)(Ts − TR ) = 0. (34) × v1 The curve 5 in Fig. 14 is drawn using Eq. (34). The curve 4 covers the area of multiplicity, and the curve 5—the area of dynamic stability. The intersection of heat production curve with the boundaries of multiplicity and limit cycles areas is noted as A and B respectively (Fig. 14). There is the area of static stability above point B. But limit cycles take place between points A and B. The farther from point A, the higher an amplitude of a limit cycle. Similarly for the case when the gas solubility increases with the temperature increasing (Q < 0) it is drawn Fig. 15 at the following values of parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 9 × 106 s−1 ; H e∞ = 6 × 104 Pa m3 /kmol; Q = −2.5 × 106 J/kmol; E = 4.3 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); = 0.18 s−1 ; P = 101000 Pa; qL = 9 × 10−4 m3 /s; = 11.1 s; L = 1300 kg/m3 ;
(1) the gas solubility is reduced as temperature increases; (2) the gas solubility increases as temperature increases. The dependence of gas solubility on temperature was described using the equation for Henry’s law constant H e = H e∞ exp(−Q/RT ), the absorption heat Q was accepted as a constant magnitude not depending on temperature if process goes at temperatures far from the boiling point of solvent. It has been shown that the gas solubility behavior has an influence on the form of heat production curve. The heat production curve has a maximum in the first case (when gas solubility is reduced as temperature increases), so the heat production is reduced at high temperatures. The heat production curve has no maximum in the second case (when gas solubility increases as temperature increases), the heat production increases at any temperature. The areas of multiplicity and trifurcation points were found. In Figs. 3,4 these areas are shown in coordinates (TR ; 1/Q 1 ) and (T ; 1/Q 1 ) for two cases of gas solubility behavior. The temperature values TR , T and parameter 1/Q 1 in trifurcation points were found. It has been obtained, that for the case of Fig. 3 the multiplicity is impossible at values 1/Q 1 < 0.3, that corresponds to values TR = 311 K, T = 348.4 K. For the case of Fig. 4, multiplicity is impossible at values 1/Q 1 < 0.5 that corresponds to values TR = 331.1 K, T = 371.7 K. The phenomenon of hysteresis in the reactor was described. The hysteresis diagram was drawn in coordinates (; T0 ) and (T; T0 ) for the case when the gas solubility is reduced as temperature increases (Fig. 5). The interval of hysteresis was found: (T0 )ign − (T0 )ext = 67.9 ◦ C. The similar diagram was drawn for the case when the gas solubility increases as temperature increases (Fig. 6). In this case the hysteresis interval is: (T0 )ign − (T0 )ext = 49.5 ◦ C. A criterion for determining whether the multiplicity is possible under the given conditions was found. The intersection points of F () with ad correspond to steady states of reactor. The function F () has two extrema which become less expressed at increasing of parameter k0 (Fig. 7) and disappear completely at certain value of k0 . Then a
4146
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
single intersection point of function F () with ad , is possible, i.e. only one steady state at any values of ad takes place at such a form of curve F (). It means that the multiplicity is impossible. By means of analysis of mathematic model the stability of steady state of the reactor was studied. The plots of gas concentration CGL and conversion dependences on dimensionless deviation from temperature TS were drawn. The types of equilibrium positions were determined as stable nodes and shown in Figs. 9 and 10 for data in Figs. 1 and 2. The conditions of static and dynamic stability of reactor were found. The reactor is dynamically unstable if it operates cyclically approaching to limit cycle. Such limit cycle is shown in Fig. 12 in coordinates (; ). The approach to limit cycle in coordinates (CGL ; t) and (; t) is presented in Fig. 13. The diagrams of reactor stability for the cases with different laws of gas solubility behavior are drawn (Figs. 14, 15). The areas of multiplicity, dynamic stability, limit cycles are shown in given diagrams.
TS
Notation
CL0 CGL ∗ CG CT L E FC He H e∞
H k k∞ kS k0 KC P qL Q Qin L Qin L r R S St T TC
liquid reactant concentration at the inlet of reactor, kmol/m3 gas dissolved concentration in liquid, kmol/m3 equilibrium gas dissolved concentration in liquid, kmol/m3 heat capacity of liquid phase, J/(kg K) activation energy, J/kmol heat transfer area, m2 Henry’s law constant, Pa m3 /kmol pre-exponential factor of Henry’s law constant, Pa m3 /kmol heat of reaction, J/kmol reaction rate constant, s−1 pre-exponential factor in the Arrhenius equation, s−1 reaction rate constant at the temperature TS , s−1 parameter k0 = k∞ exp(−1/), s−1 overall heat transfer coefficient, W/m2 K lpressure, Pa volumetric flow rate of liquid, m3 /s amount of heat, J heat entered the reactor with a liquid flow, W heat removed from the reactor with a liquid flow, W reaction rate, kmol/(m3 s) gas constant, J/kmol K upper steady-state point Stanton number, St = KCcTFLcqL L temperature,K coolant temperature, K
T0 TR Tmax
Tad VL
temperature corresponding to the upper steadystate point, K liquid temperature at the inlet of reactor, K temperature corresponding to the intersection point of heat removal line with abscissas axis +St·TC TR = T0 1+St ,K temperature corresponding to trifurcation point, K adiabatic temperature rise of a reaction mixture at full liquid reactant conversion, K volume of liquid phase, m3
Greek letters
1 , 2 ad
c ad L
overall mass-transfer coefficient between gas and liquid, s−1 0 dimensionless parameter, = RT E stoichiometric coefficients dimensionless deviation from temperature TS , = E/(RT 2S ) · (T − TS ) dimensionless adiabatic temperature rise of a reaction mixture at full liquid reactant conversion ad = Tad E/(RT 2S ) dimensionless temperature = E/(RT 20 ) · (T − T0 ) dimensionless parameter corresponding to the temperature of coolant, C = E/(RT 20 ) · (TC − T0 ) dimensionless parameter ad = Tad E/(RT 20 ) density of liquid phase, kg/m3 residence time of solution in reactor, s conversion of liquid reactant.
Abbreviations HPR HWR CSTR
heat production rate heat withdrawal rate continuous stirred tank reactor.
References Aris, R., Amundson, N.R., 1958. An analysis of chemical reactor stability and control—part I—III. Chemical Engineering Science 7, 121–155. Beskov, W.S., Charkova, T.V., Novikov, E.A., 1979. Number and stability of steady-state regimes of gas–liquid reactors. Theoretical Foundations of Chemical Technology 13, 120. Bilous, O., Amundson, N.R., 1955. Chemical reactor stability and sensitivity. A.I.Ch.E. Journal 1, 513–521. Frank-Kamenetski, D.A., 1969. Diffusion and Heat Transfer in Chemical Kinetics. Second Edition. Plenum, New York. Hoffman, L.A., Sharma, S., Luss, D., 1975. Steady-state multiplicity of adiabatic gas–liquid reactors. I. The single reactor case. A.I.Ch.E. Journal 21, 318–326. Huang, D.T.-J., Varma, A., 1981. Steady-state and dynamic behavior of fast gas–liquid reactions in non-adiabatic continuous stirred tank reactors. Chemical Engineering Journal 21, 47–57.
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147 Olsen, R.J., Epstein, I.R., 1993. Bifurcation analysis of chemical reaction mechanisms. II. Hopf bifurcation analysis. Journal of Chemical Physics 98, 2805–2822. Sharma, S., Hoffman, L.A., Luss, D., 1976. Steady-state multiplicity of adiabatic gas–liquid reactors: II. The two consecutive reactions case. A.I.Ch.E. Journal 22, 324–331. Singh, C.P.P., Shah, Y.T., 1982. The effect of gas feed temperature on the steady-state multiplicity of an adiabatic CSTR with a fast pseudofirst-order reaction. Chemical Engineering Journal 23, 101–104.
4147
Slinko, M.G., Muler, A.L., 1961. Kin. i Kat 2, 467(in Russian). van Elk, E.P., Borman, P.C., Kuipers, J.A.M., Versteeg, G.F., 1999. Modeling of gas–liquid reactors—stability and dynamic behavior of gas–liquid mass transfer accompanied by irreversible reaction. Chemical Engineering Science 54, 4869–4879. Vleeschhouwer, P.H.M., Garton, R.D., Fortuin, J.M.H., 1992. Analysis of limit cycles in an industrial oxo reactor. Chemical Engineering Science 47, 2547–2552.
Study of gas–liquid stirred reactor stability taking into account the temperature influence on gas solubility E.F. Stefoglo∗ , I.V. Kuchin, A.V. Kravtsov Institute of Coal and Coal Chemistry, Rukavishnikov, 21, 650610, Kemerovo, Russia Received 3 December 2003; received in revised form 6 April 2004; accepted 25 April 2004
Abstract The paper is devoted to the dynamic behavior and stability of gas–liquid stirred reactor taking into account the temperature influence on gas solubility. Since the rate of gas–liquid processes is very sensitive to concentration of gas reactant dissolved in liquid, even weak fluctuations of temperature can significantly influence on process pass. There are two cases of temperature influence on gas solubility are possible: (1) the solubility decreases with increasing temperature; (2) the solubility grows with increasing temperature. The first case is typical for majority of gases. The second case occurs more rarely but has a great practical importance. It takes place, for example, for the hydrogenation of many compounds in organic solvents (such as benzene, toluene, isopropyl alcohol and others). A model of gas–liquid process has been developed to demonstrate the stability of gas–liquid reactor. It has been shown that the gas solubility behavior has an influence on the form of heat production curve and therefore on the multiplicity of the steady states. The areas of multiplicity and limit cycles were found and the phenomenon of hysteresis in the reactor was shown. A criterion to determine whether the multiplicity is possible under the given conditions was found. By means of an analysis of a mathematical model the stability of steady states of the reactor was studied. 䉷 2004 Elsevier Ltd. All rights reserved. Keywords: Multiphase reactors; Dynamic simulation; Stability; Mathematical modeling; Multiplicity; Absorption
1. Introduction Multiplicity, stability and dynamic behavior of singlephase processes have been well studied in the literature. During recent decades, many papers have been devoted to this subject. The main results in this area are shown in papers (Bilous and Amundson, 1955; Aris and Amundson, 1958; Slinko and Muler, 1961; Frank-Kamenetski, 1969; Olsen and Epstein, 1993 and others) in detail. Prediction of the dynamic behavior of multiphase reactors is usually more complex and less investigated because these systems involve: (1) more than two component balances, i.e. ODEs, (2) mass transfer between the gas and
∗ Corresponding author. Tel.: 7-3842-36-55-61; fax: 7-3842-21-18-38.
E-mail address: [email protected] (E.F. Stefoglo) 0009-2509/$ - see front matter 䉷 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2004.04.041
the liquid phase and (3) gas reactant that comes into the reaction zone independently of the liquid reactant. The existing investigations (Hoffman et al., 1975; Sharma et al., 1976; Beskov et al., 1979; Huang and Varma, 1981; Singh and Shah, 1982; Vleeschhouwer et al., 1992; van Elk et al., 1999) do not take into account that the changing temperature leads to the changes in the gas solubility in the liquid phase. However, very often the real processes are carried out at the unsteady temperature. Changes in the temperature lead to the changes in gas solubility in liquid phase. Dependence of gas solubility on the temperature really can be very significant. The temperature change in 10 ◦ only can lead sometimes to the change of gas solubility several times (such systems as hydrogen–acetone, hydrogen–benzene, etc.). Since the concentration of gas reactant dissolved in liquid is, as a rule, much lower than the liquid reactant concentration, the gas concentration can influence the reaction rate to a great
4138
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
extent. Thus, in the study of thermal stability of a reactor it is very important to take into account the effect of temperature on gas reactant content in the liquid. The dynamics of gas–liquid process are considered in this paper by means of a model accounting for the temperature influence on gas solubility. First-order reaction (r = kC GL ) is considered of the type G(gas) + L(liquid reactant) → P (product) carried out in continuous stirred tank reactor (CSTR).
2. Mathematical model of CSTR taking into account temperature influence on gas solubility For steady-state operating gas–liquid reactor the material and heat balance are: for a gas reactant: ∗ VL (CG − CGL ) − VL 1 kC GL − qL CGL = 0,
(1)
Two cases of temperature influence on gas solubility are possible: 1. solubility decreases as temperature increases. 2. solubility increases as temperature increases. The first case is typical for majority of gases. The Henry’s law constant depends on temperature according to equation H e = H e∞ exp(−Q/RT ). The second case is not so common but is of great practical importance. It takes place, for example, with hydrogen dissolution in organic solvents (benzene, toluene, isopropyl alcohol and others). In this case, the Henry’s law constant must decrease with increasing temperature. This is possible if the heat of absorption Q in expression H e = H e∞ exp(−Q/RT ) is formally considered as a negative magnitude. In the quite wide range of temperature far from the boiling point of solvent the absorption heat Q may be accepted as a constant magnitude not depending on temperature. Introducing the Stanton number as
for a liquid reactant: qL CL0 − qL CL − VL 2 kC GL = 0,
(2)
(3)
where Qin L = qL L T0 CT L is the heat entered the reactor by the liquid flow, [W], Qout L = qL L CT L T the heat removal from the reactor with a liquid flow, [W], QR = VL H r is the heat production as a result of chemical reaction, [W], Kc Fc (T − Tc ) the heat removal from the reactor through heat exchange surface, [W]. ∗ = P /H e), From Eqs. (1)–(2) and Henry’s law (CG Eq. (3) takes the form: Kc Fc Tad · = (T − T0 ) + (T − Tc ). L CT L qL
Kc Fc , L C T L q L
we write the right-hand side of Eq. (4) as
heat balance equation: out Qin L − QL + QR − Kc Fc (T − Tc ) = 0,
St =
(4)
(1 + St)(T − TR ), where TR =
T0 + St · Tc . 1 + St
(5)
3. Results of model solution. Multiplicity of steady states Heat production curve Q1 and heat removal line Q2 are shown in Fig. 1 for the case when gas solubility is
The left-hand side of Eq. (4) represents the heat production (Q1 ) as a result of chemical reaction, while the righthand side represents the heat removal (Q2 ) with outgoing liquid flow and heat exchange. In Eq. (4) the following magnitudes are used:
Tad =
H C L0 — adiabatic heating, L C T L
k∞ e−E/RT · P 2 CL0 H e∞ e−Q/RT ( + k∞ e−E/RT + 1) — conversion of liquid reactant,
=
is the residence time of the solution in the reactor, [s], H e =H e∞ exp(−Q/RT ) is the dependence of Henry’s law constant on the temperature; Q is the heat of absorption, [J/kmol].
Fig. 1. Heat production curves and heat removal line, steady-state multiplicity in CSTR for a the first-order reaction. The gas solubility is reduced as temperature increases.
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
Fig. 2. Heat production curve and heat removal line for the first-order reaction in CSTR. The gas solubility increases as temperature increases.
reduced with temperature increasing. The plot is drawn at the following values of process parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 1 × 106 s−1 ; H e∞ = 1 × 107 Pa m3 /kmol; Q=3×106 J/kmol; E=4.9×107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); = 0.07 s−1 ; P = 101000 Pa; qL = 8 × 10−6 m3 /s; = 1.25 × 103 s; L = 1300 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.
Fig. 3. Areas of steady states multiplicity in gas–liquid reactor in terms of temperatures T and TR = (Q 1 · T − Q1 )/Q . First-order reaction. The gas solubility is reduced as temperature increases. Calculations were made on data of Fig. 1.
At the inflection point (TF in Figs. 1 and 2) a slope of tangent to heat production curve is maximum. If the heat removal line has a slope greater than the maximum slope, then only one intersection point of heat production curve with heat removal line is possible. The multiplicity is not possible in such conditions. The tangency point TF with the maximum slope of tangent is a trifurcation point. Using expression (4) we find the slope of tangent to heat production curve at trifurcation point:
E−Q · P 2 k∞ e−E/RT max RT ( + k∞ e−E/RT max + 1) − k∞ e−E/RT max RTE2 2 1 + St d(Q1 )max max max = , = Tad dT H e∞ e−Q/RT max CL0 ( + k∞ e−E/RT max + 1)2 It has been seen from Fig. 1, that at T > 430 K the heat production is reduced because the gas solubility decreases with temperature increasing. Two positions of heat removal line are shown in Fig. 1a and b—the tangents to heat production curve. Both lines have the same slope and intersect an abscissa axis at points TRa and TRb . Three areas can be indicated in Fig. 1: (1) TR < TRa : in this area only one stable steady state exists with low value of liquid reactant conversion. (2) TRa < TR < TRb : in this area three steady-states occur—two stable states with low and high conversion, one unstable state with middle conversion. (3) TR > TRb : in this area only one stable steady state with high conversion exists. Heat production curve Q1 and heat removal line Q2 are shown in Fig. 2 for the case when with temperature increasing the gas solubility increases (Q < 0). The plot is drawn at the following values of parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 2.2 × 106 s−1 ; H e∞ = 3.5 × 105 Pa m3 /kmol; Q = −2.5 × 106 J/kmol; E = 5 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); = 0.2 s−1 ; P = 101000 Pa; qL = 2 × 10−4 m3 /s; = 50 s; L = 1300 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 275 K; KC = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.
4139
(6)
where Tmax is the temperature corresponding to trifurcation point TF. To find the temperature Tmax corresponding to trifurcation point, it is necessary to solve numerically the following equation: d2 (Q1 ) = 0. (7) dT 2 For data of Fig. 1 Tmax is equal to 348.4 K. Thus, a criterion of a single stable steady-state existence is 1 + St d(Q1 )max > , Tad dT
(8)
1 )max is calculated using Eq. (6). where d(QdT It should be noted, that value 1/Q 1 = 1/(dQ1 /dT ) is the maximum increase in reactor temperature if full conversion of the liquid reactant takes place under the non-adiabatic conditions. The areas where multiplicity can take place are shown in Fig. 3 in coordinates (TR ; 1/Q 1 ) and (T; 1/Q 1 ). As can be seen from Fig. 3, a single, stable steadystate takes place at 1/Q 1 < 0.3 that corresponds to values TR = 311 K, T = Tmax = 348.4 K. From Fig. 3, it follows that at a value 1/Q 1 = 0.7 the multiplicity takes place within the temperature intervals 227.8 < TR < 291.2 K and 310 < T < 387.7 K.
4140
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
Fig. 5. The hysteresis in CSTR for the case when gas solubility decreases as temperature increases (Q > 0). Calculations were made on data of Fig. 1. Fig. 4. Areas of multiplicity in gas–liquid reactor in terms of temperatures T and TR = (Q 1 · T − Q1 )/Q . First-order reaction. The gas solubility increases as temperature increases. Calculations were made on data of Fig 2.
Fig. 4 is drawn similarly to Fig. 3 but for the case in which the gas solubility increases with increasing temperature. We have one stable steady state in this case at 1/Q 1 < 0.5 that corresponds to values TR = 331.1 K, T = Tmax = 371.7 K.
4. Influence of gas solubility on the hysteresis in CSTR Let us consider the reactor operating in a stable steadystate with low value of liquid reactant conversion (Fig. 1). As the temperature T0 of the liquid flowing to the reactor inlet, the heat removal line moves in parallel to itself to the right and at the certain moment will coincide with a tangent to heat production curve (line b in Fig. 1). As T0 increases further, only one steady-state characterized by the high conversion will take place. Thus, at a given value of T0 called the ignition temperature (T0 )ign , a strong increase in the reaction rate has occurred and the steady state has moved to the upper position. For the opposite case, i.e., for a reduction in the temperature T0 of liquid flowing to the reactor inlet, the process occurs at a high conversion until the reaction is extinguished at a certain temperature (T0 )ext . The hysteresis phenomena means that the ignition and extinction temperatures are dissimilar: (T0 )ign = (T0 )ext . The hysteresis diagram (Fig. 5) is calculated from the data of Fig. 1. The hysteresis interval in this case is equal: (T0 )ign − (T0 )ext = 67.9 ◦ C. A similar diagram is obtained for the case in which the gas solubility increases with temperature (Q < 0). The results of this calculation are presented in Fig. 6. Calculations were made on the data of Fig. 2. In this case the hysteresis interval is constant: (T0 )ign − (T0 )ext = 49.5 ◦ C.
Fig. 6. The hysteresis in CSTR for case when with gas solubility increases as temperature increases (Q > 0). Calculations were made on data of Fig. 2.
5. Criterion of steady-states multiplicity in the reactor We now develop a criterion indicating the possibility of multiplicity in gas–liquid reactor, using the approach proposed by Slinko (1961). From Eq. (4) we have heat production (HPR): Tad · , heat removal (HWR): (T − T0 ) +
Kc F c (T − Tc ). L CT L qL
In the steady-state HWR = HPR or HWR/HPR = 1. Let us introduce the following dimensionless values:
= E/(RT 20 )(T − T0 );
T =
RT 20 + T0 ; E
RT 0 = ; E L ; exp −
L = Q/E; S = P /(H e0 CL0 ); H e0 = H e∞ k = k0 exp ; k0 = k∞ exp(−1/). 1 +
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
4141
The differentiation of function F () gives dF d A B + (−q)2 k0 exp 1+ − B(1 − L) (1+) , = (1−L) 2 S k0 exp 1+
F () ≡
(11) where A = (1 + U k 0 · ), B = + k0 exp Fig. 7. The function F () at various values of parameter k0 . Calculations were made on data of Fig. 1.
Then Eq. (4) can be rewritten in the form: E HWR Tad 2 HPR RT 0 Kc Fc E ) + (T − T ) (T − T 0 c 2 C q RT L TL L . = 0
≡F (),
(13)
k0 exp
− B(1 − L) 1 + (14)
By substituting Eq. (13) in Eq. (14), after some transformations we discover: k0 exp 1 + ( − q) =− + k exp + 1 0 1 + (1 + )2
×L − − 1 − − 1.
(15)
(10) U ∗ k0 c 1+U ∗ k0 ;
c = E/(RT 20 ) · (Tc − T0 ); U ∗ = Kc Fc 2 k0 L CT L qL ; ad = Tad E/(RT 0 ). At steady states, the relation HWR/HPR = 1. The intersection points of curve F () and line ad correspond to the steady states of the reactor. We examine the function F () to define when the multiplicity is possible. In Fig. 7 a plot of the function F () at various values of parameter k0 is shown. For different values of k0 there can exist three solutions or just one, depending on the value ad . In Fig. 7 ad = 20.7, and only at k0 = 2.6 three intersection points of function F () with ad are possible. At other values of parameter k0 only one intersection point exists. With an increase in the parameter k0 the extrema of function F () become less expressed and disappear completely at certain k0 value. For such a form of the curve F () only one intersection point of function F () and ad can occur. This means that there is a single steady state for any values of ad , and thus multiplicity is impossible. Therefore, it is necessary to formulate a condition when the function F () loses the extrema. where q =
( − q)
(1 + )2 =0.
If the left-hand side of Eq. (9) note as ad and the above dimensionless complexes introduce into the right-hand side, Eq. (9) takes the form: ( − q)(1 + U ∗ k0 ) · + k0 exp 1+ + 1 ad = (1−L) k0 exp 1+ 2 · S
+ 1. 1 +
To find the extrema of function F () we equate F () to zero. Taking into account that A = 0 and the denominator of expression (11) is not zero, we find: B+
(9)
(12)
The left-hand side of expression (15) is some function of (we note it as f1 ()). The right-hand side of this expression we note as f2 (). Then expression (15) takes the following form: f1 () = f2 (). If f1 () and f2 () intersect each other at some points 1 and 2 , then derivative F is equal to zero at these points, therefore the function F has extrema. If f1 () = f2 (), the functions f1 and f2 do not intersect each other, F = 0 and function F has no extrema, therefore the multiplicity is impossible. Thus, the problem of determination of function F () extrema reduces to the finding of intersection points of the curves f1 () and f2 (). To define whether the intersection of functions f1 and f2 takes place, the following method can be used. We draw two parallel tangents to functions f1 and f2 by such way that the line passing through tangency points 1 and 2 to be perpendicular to the tangents (Fig. 8). Values 1 and 2
4142
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
characterized by the high value of conversion. We note this steady state by point S. Using the expansion proposed by Frank–Kamenetski (1969) we find:
Fig. 8. Functions f1 and f2 dependences on .
corresponding to points 1 and 2 can be found by numerical solution of the following equations set: f1 (1 ) = f2 (2 ), f2 (2 ) − f1 (1 ) = −1. f1 (1 ) 2 − 1
(16)
The first equation means that the tangents to functions f1 and f2 at points 1 and 2 are parallel (their derivatives are equal). The second equation means that the line passing through tangency points 1 and 2 is perpendicular to tangents (if lines are perpendicular, the product of their angular factors is equal to −1). Calculation of 1 and 2 allows to we find values of functions f1 (1 ) and f2 (2 ). The condition when the functions f1 and f2 are not intersected, and therefore the multiplicity is impossible is f1 (1 ) > f2 (2 )
or
2 > 1 .
(17)
6. Stability of steady states of reactor In practice, we deal with non-stationary operation of continuous reactor because the various deviations of process parameters are possible. Therefore it is necessary to know when these deviations do not lead to the change of process mode, i.e. to know the conditions of stable reactor operation. Under the non-stationary conditions the mass balance equation for gas reactant is: dCL ∗ − CGL ) − 1 k CGL − CGL , = (CG dt
(18)
where t = ϑ/ is the dimensionless time; ϑ the current time, ∗ = P /H e the equilibrium [s]; the residence time, [s]; CG concentration of gas dissolved in liquid. On the basis of expression (4) we write an equation for temperature in reactor: dT = (T0 − T ) − St(T − Tc ) + Tad · . dt
(19)
As described above (Figs. 1 and 2), there are no more than three intersection points of heat production curve and heat removal line. One point corresponds to an unstable steady state and two points—to stable ones. The upper stable steady-state represents the special interest because it is
k = kS exp(),
(20)
H e = H eS exp( · L),
(21)
where H es is the Henry’s law constant at the temperature, Ts of the upper stable steady state, [Pa m3 /kmol]; the dimensionless deviation from temperature Ts , = E/(RT 2s )(T − Ts ); L = Q/E the dimensionless parameter; Q the heat of absorption, [J/kmol]; E the activation energy, [J/kmol]; ks the reaction rate constant at the temperature Ts , [s−1 ]; By dividing of Eq. (18) on Eq. (19), using dependences (20), (21) and introducing dimensionless parameters, we find H e Pe·L −CGL −CGL (v1 ks e +1) dCGL s = , d −(1 + St)+ad v2 ksve1 CLCGL −ad s
(22)
0
where ad = Tad E/(RT 2S ) is the dimensionless parameter, s the conversion of liquid reactant at the temperature Ts ; St = Kc Fc /(L CT L qL ) the Stanton’s criterion; 1 , 2 the stoichiometric coefficients. From Eq. (22) the dependence of gas concentration CGL on the dimensionless temperature deviation from a steady state was found. Knowing the dependence CGL on , it is easy to find the function on with the aid of relation:
= ks e · CGL /CL0 .
(23)
In Fig. 9 the dependence on in combination with heat production and heat removal lines are shown. If as a result of some disturbance, the process deviates from stable state (point S) to some state (point A), then the process will return to S by the trajectories presented in Fig. 9. A similar plot for the case when the solubility of gas increases with temperature increasing (Q < 0) is shown in Fig. 10. Solution of Eq. (22) is the integral curves of the following system: dCGL P = − CGL − CGL (v1 ks e + 1) dt H eS e·L ≡1 (CGL , ), v2 ks e CGL d =−(1 + St) + ad − ad s dt v1 CL0 ≡2 (CGL , ).
(24)
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
4143
Fig. 11. Areas corresponding to various types of equilibrium position.
= 0 in the upper steady state, the coefficients of the linearized equations willbe equal: Fig. 9. The dependence and CGL on for the first-order reaction in gas–liquid CSTR. The solubility of gas increases as temperature decreases (Q > 0). Calculations were made on data of Fig. 1.
a11 =
*1 = − − v1 ks − 1, *CGL
a12 =
*1 P =− L − v1 ks · CGLs , H es *
a21 =
*2 ad v2 ks = , v1 CL0 *CGL
a22 =
*2 ad v2 ks · CGLs = −(1 + St) + . v1 CL0 *
The type of equilibrium position is determined by values of known criterions: = −(a11 + a22 ), = a11 a22 − a12 a21 . (26)
Fig. 10. The dependence and CGL on for the first-order reaction in gas–liquid CSTR. The solubility of gas increases as temperature increases (Q < 0). Calculations were made on data of Fig. 2.
The points of phase plane corresponding to equilibrium position (steady state of the reactor) satisfy to the conditions: dCGL d = = 0. dt dt Then
1 (CGLS , S ) = 2 (CGLS , S ) = 0,
(25)
where CGLS , S are values of gas concentration and temperature deviation under equilibrium condition (at the temperature TS ). To study the process stability it is convenient to linearize Eqs. (24), near the temperature TS . Taking into account that
The positive values of parameters and correspond to stable equilibrium positions (Fig. 11). Then the conditions of stability in terms of expressions for coefficients a11 , a12 , a21 , a22 , are: ad v2 ks · CGLs (− − v1 ks − 1) · −(1 + St) + v1 CL0 P ad v2 ks − L − v1 ks · CGLs · > 0, (27) H es v1 CL0 ad v2 ks · CGLs ( + v1 ks + 1) − −(1 + St) + v1 CL0 >0 (28) or using expression (23) and taking into account that S =0:
E v 2 S (− − v1 ks − 1) · −(1 + St) + RT 2S v1 P E Tad v2 kS − L − v1 S CL0 · > 0, (27a) H es RT 2S v1 CL0
4144
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
Fig. 12. Limit cycle, dependence on for the first-order reaction in gas–liquid CSTR.
E v2 S ( + v1 ks + 1) − −(1 + St) + RT 2S v1
> 0. (28a)
The areas corresponding to various types of equilibrium positions in coordinates , are shown in Fig. 11. The point A corresponds to values and calculated on data of Fig. 1, the point B—on data of Fig. 2. Eq. (27) represents a static stability condition, and Eq. (28), the dynamic stability condition. If conditions (27) and (28) are met, the process is stable in steady state. If condition (27) is met and condition (28) is not met, the reactor is dynamically unstable and operates cyclically in the limit cycle around the point S corresponding to the upper steady state (Fig. 12). The plot is drawn at following values of parameters: VL =0.01 m3 ; H =1.7×108 J/kmol; k∞ =1×106 s−1 ; H e∞ = 2.5 × 107 Pa m3 /kmol; Q = 2 × 107 J/kmol; E = 4.9 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); =0.05 s−1 ; P =101000 Pa; qL =9×10−4 m3 /s; =11.1 s; L = 1200 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.The transition to a limit cycle will take place as a result of process deviation from a steady state to point A. The change of process variables (gas dissolved concentration, conversion of liquid reactant, temperature) in the reactor has a cyclical character with a certain frequency and amplitude. The oscillations of values CGL and corresponding to the case of Fig. 12 are shown as a function of time in Fig. 13. We draw the diagram of stability of the given reactor. From expression (4) we find a relation between S and TS :
S =
1 + St (TS − TR ), Tad
where S =
k∞ e−E/RT S ·P 2 −Q/RT S (+k∞ e−E/RT S +1) CL0 H e∞ e
(29)
Fig. 14. Stability diagram for CSTR. Gas solubility increases as temperature increases (Q < 0). (1) zero conversion; (2) full conversion; (3) heat production; (4) area of multiplicity; (5) area of limit cycles.
Dependence (30) is shown as a curve 3 in Fig. 14. Calculations were made at following values of parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 9 × 106 s−1 ; H e∞ = 5 × 105 Pa m3 /kmol; Q = 9 × 106 J/kmol; E = 4.9 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); =0.05 s−1 ; P =101000 Pa; qL =9×10−4 m3 /s; =11.1 s; L = 1200 kg/m3 ; CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K. At zero conversion (S = 0) we have from Eq. (30): TS − TR = 0.
(30)
(31)
This condition corresponds to line 1 in Fig. 14. Line 2 corresponds to full conversion S = 1: TS − TR =
.
Rewriting Eq. (29), we find:
Tad TS − TR = . 1 + St S
Fig. 13. Oscillations of values CGL and during time.
Tad . 1 + St
(32)
Eq. (27a) is a condition of static stability. To find the boundary of static stability region, we equate the lefthand side of Eq. (27a) to zero, substituting instead of S
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
4145
CT L = 1 × 103 J/(kg K); T0 = 273 K; Kc = 95 W/(m2 K); Fc = 0.081 m2 ; Tc = 253 K.
7. Conclusions Since dissolved gas concentration is a magnitude of several orders lower than liquid reactant concentration, the rate of gas–liquid processes is very sensitive to concentration of gas reactant dissolved in liquid, and so even weak fluctuations of temperature can significantly influence on process pass. In the given paper the dynamics of gas–liquid process in the CSTR with taking into account the temperature influence on gas solubility has been considered for two cases: Fig. 15. Stability diagram for CSTR. Gas solubility reduces as temperature increases (Q > 0). (1) zero conversion; (2) full conversion; (3) heat production; (4) area of multiplicity; (5) area of limit cycles.
the right-hand side of Eq. (29): (− − v1 ks − 1) E v2 (1 + St)(Ts − TR ) × −(1 + St) + v1 RT 2s P 1 + St L − CL0 (Ts − TR ) − − H es Tad E Tad v2 ks × = 0. RT 2s v1 CL0
(33)
Using Eq. (33) we draw a curve 4 in Fig.14, which covers the statically unstable area of reactor operation (area of multiplicity). As a result of similar operations with condition (28a) of dynamic stability, we obtain E ( + v1 ks + 1) − −(1 + St) + RT 2s v2 (1 + St)(Ts − TR ) = 0. (34) × v1 The curve 5 in Fig. 14 is drawn using Eq. (34). The curve 4 covers the area of multiplicity, and the curve 5—the area of dynamic stability. The intersection of heat production curve with the boundaries of multiplicity and limit cycles areas is noted as A and B respectively (Fig. 14). There is the area of static stability above point B. But limit cycles take place between points A and B. The farther from point A, the higher an amplitude of a limit cycle. Similarly for the case when the gas solubility increases with the temperature increasing (Q < 0) it is drawn Fig. 15 at the following values of parameters: VL = 0.01 m3 ; H = 1.7 × 108 J/kmol; k∞ = 9 × 106 s−1 ; H e∞ = 6 × 104 Pa m3 /kmol; Q = −2.5 × 106 J/kmol; E = 4.3 × 107 J/kmol; CL0 = 2 kmol/m3 ; R = 8314 J/(kmol K); = 0.18 s−1 ; P = 101000 Pa; qL = 9 × 10−4 m3 /s; = 11.1 s; L = 1300 kg/m3 ;
(1) the gas solubility is reduced as temperature increases; (2) the gas solubility increases as temperature increases. The dependence of gas solubility on temperature was described using the equation for Henry’s law constant H e = H e∞ exp(−Q/RT ), the absorption heat Q was accepted as a constant magnitude not depending on temperature if process goes at temperatures far from the boiling point of solvent. It has been shown that the gas solubility behavior has an influence on the form of heat production curve. The heat production curve has a maximum in the first case (when gas solubility is reduced as temperature increases), so the heat production is reduced at high temperatures. The heat production curve has no maximum in the second case (when gas solubility increases as temperature increases), the heat production increases at any temperature. The areas of multiplicity and trifurcation points were found. In Figs. 3,4 these areas are shown in coordinates (TR ; 1/Q 1 ) and (T ; 1/Q 1 ) for two cases of gas solubility behavior. The temperature values TR , T and parameter 1/Q 1 in trifurcation points were found. It has been obtained, that for the case of Fig. 3 the multiplicity is impossible at values 1/Q 1 < 0.3, that corresponds to values TR = 311 K, T = 348.4 K. For the case of Fig. 4, multiplicity is impossible at values 1/Q 1 < 0.5 that corresponds to values TR = 331.1 K, T = 371.7 K. The phenomenon of hysteresis in the reactor was described. The hysteresis diagram was drawn in coordinates (; T0 ) and (T; T0 ) for the case when the gas solubility is reduced as temperature increases (Fig. 5). The interval of hysteresis was found: (T0 )ign − (T0 )ext = 67.9 ◦ C. The similar diagram was drawn for the case when the gas solubility increases as temperature increases (Fig. 6). In this case the hysteresis interval is: (T0 )ign − (T0 )ext = 49.5 ◦ C. A criterion for determining whether the multiplicity is possible under the given conditions was found. The intersection points of F () with ad correspond to steady states of reactor. The function F () has two extrema which become less expressed at increasing of parameter k0 (Fig. 7) and disappear completely at certain value of k0 . Then a
4146
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147
single intersection point of function F () with ad , is possible, i.e. only one steady state at any values of ad takes place at such a form of curve F (). It means that the multiplicity is impossible. By means of analysis of mathematic model the stability of steady state of the reactor was studied. The plots of gas concentration CGL and conversion dependences on dimensionless deviation from temperature TS were drawn. The types of equilibrium positions were determined as stable nodes and shown in Figs. 9 and 10 for data in Figs. 1 and 2. The conditions of static and dynamic stability of reactor were found. The reactor is dynamically unstable if it operates cyclically approaching to limit cycle. Such limit cycle is shown in Fig. 12 in coordinates (; ). The approach to limit cycle in coordinates (CGL ; t) and (; t) is presented in Fig. 13. The diagrams of reactor stability for the cases with different laws of gas solubility behavior are drawn (Figs. 14, 15). The areas of multiplicity, dynamic stability, limit cycles are shown in given diagrams.
TS
Notation
CL0 CGL ∗ CG CT L E FC He H e∞
H k k∞ kS k0 KC P qL Q Qin L Qin L r R S St T TC
liquid reactant concentration at the inlet of reactor, kmol/m3 gas dissolved concentration in liquid, kmol/m3 equilibrium gas dissolved concentration in liquid, kmol/m3 heat capacity of liquid phase, J/(kg K) activation energy, J/kmol heat transfer area, m2 Henry’s law constant, Pa m3 /kmol pre-exponential factor of Henry’s law constant, Pa m3 /kmol heat of reaction, J/kmol reaction rate constant, s−1 pre-exponential factor in the Arrhenius equation, s−1 reaction rate constant at the temperature TS , s−1 parameter k0 = k∞ exp(−1/), s−1 overall heat transfer coefficient, W/m2 K lpressure, Pa volumetric flow rate of liquid, m3 /s amount of heat, J heat entered the reactor with a liquid flow, W heat removed from the reactor with a liquid flow, W reaction rate, kmol/(m3 s) gas constant, J/kmol K upper steady-state point Stanton number, St = KCcTFLcqL L temperature,K coolant temperature, K
T0 TR Tmax
Tad VL
temperature corresponding to the upper steadystate point, K liquid temperature at the inlet of reactor, K temperature corresponding to the intersection point of heat removal line with abscissas axis +St·TC TR = T0 1+St ,K temperature corresponding to trifurcation point, K adiabatic temperature rise of a reaction mixture at full liquid reactant conversion, K volume of liquid phase, m3
Greek letters
1 , 2 ad
c ad L
overall mass-transfer coefficient between gas and liquid, s−1 0 dimensionless parameter, = RT E stoichiometric coefficients dimensionless deviation from temperature TS , = E/(RT 2S ) · (T − TS ) dimensionless adiabatic temperature rise of a reaction mixture at full liquid reactant conversion ad = Tad E/(RT 2S ) dimensionless temperature = E/(RT 20 ) · (T − T0 ) dimensionless parameter corresponding to the temperature of coolant, C = E/(RT 20 ) · (TC − T0 ) dimensionless parameter ad = Tad E/(RT 20 ) density of liquid phase, kg/m3 residence time of solution in reactor, s conversion of liquid reactant.
Abbreviations HPR HWR CSTR
heat production rate heat withdrawal rate continuous stirred tank reactor.
References Aris, R., Amundson, N.R., 1958. An analysis of chemical reactor stability and control—part I—III. Chemical Engineering Science 7, 121–155. Beskov, W.S., Charkova, T.V., Novikov, E.A., 1979. Number and stability of steady-state regimes of gas–liquid reactors. Theoretical Foundations of Chemical Technology 13, 120. Bilous, O., Amundson, N.R., 1955. Chemical reactor stability and sensitivity. A.I.Ch.E. Journal 1, 513–521. Frank-Kamenetski, D.A., 1969. Diffusion and Heat Transfer in Chemical Kinetics. Second Edition. Plenum, New York. Hoffman, L.A., Sharma, S., Luss, D., 1975. Steady-state multiplicity of adiabatic gas–liquid reactors. I. The single reactor case. A.I.Ch.E. Journal 21, 318–326. Huang, D.T.-J., Varma, A., 1981. Steady-state and dynamic behavior of fast gas–liquid reactions in non-adiabatic continuous stirred tank reactors. Chemical Engineering Journal 21, 47–57.
E.F. Stefoglo et al. / Chemical Engineering Science 59 (2004) 4137 – 4147 Olsen, R.J., Epstein, I.R., 1993. Bifurcation analysis of chemical reaction mechanisms. II. Hopf bifurcation analysis. Journal of Chemical Physics 98, 2805–2822. Sharma, S., Hoffman, L.A., Luss, D., 1976. Steady-state multiplicity of adiabatic gas–liquid reactors: II. The two consecutive reactions case. A.I.Ch.E. Journal 22, 324–331. Singh, C.P.P., Shah, Y.T., 1982. The effect of gas feed temperature on the steady-state multiplicity of an adiabatic CSTR with a fast pseudofirst-order reaction. Chemical Engineering Journal 23, 101–104.
4147
Slinko, M.G., Muler, A.L., 1961. Kin. i Kat 2, 467(in Russian). van Elk, E.P., Borman, P.C., Kuipers, J.A.M., Versteeg, G.F., 1999. Modeling of gas–liquid reactors—stability and dynamic behavior of gas–liquid mass transfer accompanied by irreversible reaction. Chemical Engineering Science 54, 4869–4879. Vleeschhouwer, P.H.M., Garton, R.D., Fortuin, J.M.H., 1992. Analysis of limit cycles in an industrial oxo reactor. Chemical Engineering Science 47, 2547–2552.