Ind. Eng. Chem. Res. 2007, 46, 1197-1207
1197
Nonlinear Behavior of Reactor-Separator Networks: Influence of the Energy Balance Formulation K.-P. Zeyer,*,† A. A. Kulkarni,† A. Kienle,†,‡ M. Vasudeva Kumar,§ and S. Pushpavanam§ Max-Planck-Institut fu¨r Dynamik komplexer technischer Systeme, Sandtorstr. 1, D-39106 Magdeburg, Germany, Lehrstuhl fu¨r Automatisierungstechnik/Modellbildung, Otto-Von-Guericke-UniVersita¨t, UniVersita¨tsplatz 2, D-39106 Magdeburg, Germany, and Department of Chemical Engineering, Indian Institute of Technology, Chennai (Madras) 600036, India
The dynamic behavior of a stand alone flash with a constant heating rate and a corresponding reactorseparator system with recycle is analyzed by means of different model formulations which are frequently applied. It is shown that the predictions of the different model formulations are rather close for the stand alone flash. In particular, in all cases, a stable steady state is predicted. In contrast to this, it is shown that the dynamics of the corresponding reactor-separator system can depend on the specific model formulation. For a rigorous model with a dynamic energy balance, a new type of bifurcation is found, where the eigenvalues have a pole with a change of sign. Similar effects are reported for systems with negligible and finite transportation delay of the recycle. For the latter, multiple periodic solutions are found for a model with a quasistatic energy balance, which will disappear because of the new type of bifurcation if a rigorous dynamic energy balance is applied. 1. Introduction Modeling of evaporative separators with mass and energy balances, only, often leads to differential algebraic systems (DAE) with differential index of 2. To overcome the problems which are associated with the numerical solution of index 2 DAEs either a formal index reduction can be applied or the model can be reformulated by modifying the underlying assumptions. Typical modifications leading to index 1 problems are either the assumption of a quasistatic energy balance or the inclusion of some simple fluid dynamics. In the present paper, the dynamic behavior predicted by these different model formulations are analyzed and compared with each other. First, focus is on an isolated single stage flash process with given heat input. The well-known adiabatic flash is included as a special case. It is shown that the dynamic behavior of the different model formulations for a stand alone flash are rather close. In particular, in all cases, a stable steady state is predicted. Second, a simple reactor-separator system with recycle is considered. The reactor is a continuous stirred tank reactor (CSTR), whereas the separator is again a constantly heated flash. In a previous paper, we have shown that the recycle in such a system can induce instability and multiplicity of steady states.1 In the present paper, it is shown that, although the stand alone flash is rather insensitive to the specific model formulation, the dynamics of the coupled system can under certain conditions to be specified in this paper depend on the specific model formulation applied. In particular, a new type of bifurcation is found where the eigenvalues have a pole with a change of sign. The analysis is based on a model with negligible transportation delay of the recycle. Afterward, results are also extended for a finite transportation delay of the recycle. * To whom correspondence should be addressed. E-mail: zeyer@ mpi-magdeburg.mpg.de. Fax: +49/(0)391/6110-516. Phone: +49/(0)391/6110-187. † Max-Planck-Institut fu¨r Dynamik komplexer technischer Systeme. ‡ Otto-von-Guericke-Universita¨t. § Indian Institute of Technology.
Figure 1. Schematic diagrams. (a) Stand alone flash. (b) Reactor-separator with liquid recycle.
The present paper continues our fundamental research on the nonlinear behavior of reactor-separator systems. The main emphasis in our previous work was on different reaction systems such as a first-order isothermal reaction,1 a first-order exothermic reaction,2 and an autocatalytic isothermal reaction.3 In all three cases, negligible delay in the recycle loop and ideal thermodynamics were assumed. More recently, also the influence of delay4,5 and the influence of nonideal thermodynamics for azeotropic mixtures6,7 were studied. In all cases, a simple single stage flash was considered as a separator. Further, isothermal isobaric operation of the flash was assumed in most cases and thus the energy balance of the flash was neglected. In contrast to this, in refs 1 and 7, also a nonisothermal flash was considered. Consequently, the energy balance had to be taken into account. In a first step, a quasistatic energy balance was
10.1021/ie0607595 CCC: $37.00 © 2007 American Chemical Society Published on Web 12/24/2006
1198
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007
assumed. In the present paper, the influence of a more detailed energy balance formulation is studied in detail. 2. Stand Alone Flash 2.1. Model Equations. First, the focus is on a simple single stage flash process with a binary nonreactive mixture as illustrated in Figure 1a. The flash is operated at constant pressure pfl, constant heat input Qfl, and constant liquid holdup Mfl. Thermodynamic equilibrium between the vapor and the liquid phase is assumed. With these assumptions, the material balance of component A reads
Mfl
dxe ) Fz - Vye - Lxe dt
(1)
where z , xe, and ye are the mole fractions of A in the feed and the liquid and vapor phases, respectively. V and L are the vapor and liquid outflows from the flash. Since the molar holdup of the flash Mfl is constant, the total material balance is
0)F-V-L
(2)
The energy balance of the flash is obtained as
dH ) Fcp(TR - Tfl) - V∆h + Qfl dt
(4)
Here, TR is the temperature of the feed to the flash. In the remainder, a linear dependency of ∆h the enthalpy of vaporization of the mixture on the vapor-phase concentration ye is assumed, i.e.
∆h ) ye∆hA + (1 - ye)∆hB
(5)
where ∆hA and ∆hB are the enthalpies of vaporization of the pure components A and B, which are assumed to be constant. The enthalpy H on the left-hand side of the energy balance (eq 3) is the total enthalpy of the flash. For moderate pressures, the molar holdup of the vapor phase is negligible compared to the liquid phase, and the total enthalpy of the flash H can be approximated by the total enthalpy of the liquid phase, only. In general, H is a function of temperature, pressure, and composition. The time derivative of H in eqs 3 and 4 thus follows from the total differential
dt
∂H dTfl ∂Tfl dt
+
∂H dpfl ∂pfl dt
2
+
∂H dni
∑ i)1 ∂n
where cp is the total molar heat capacity of the liquid phase of the flash and hl,i denotes the molar enthalpy of component i of the liquid phase of the flash. For the analytical treatment below, for simplicity, hl,A ) hl,B is further assumed leading to
Mflcp
dTfl ) Fcp(TR - Tfl) - V∆h + Qfl dt
(8)
The vapor-liquid equilibrium is described by
ye ) KA(Tfl, xe)xe
(9)
1 ) KA(Tfl, xe)xe + KB(Tfl, xe)(1 - xe)
(10)
Assuming ideal liquid and vapor phases, the equilibrium constants follow from
pis(Tfl) , i ) A, B pfl
(11)
(3)
where hf, hv, and hl are the total molar enthalpies of the feed of the flash and the vapor and the liquid outlets, respectively. In view of hv ) hl + ∆h, hf - hl ) cp(TR - Tfl) and the total material balance (eq 2), the energy balance (eq 3) reduces to
)
dxe dTfl + Mfl(hl,A - hl,B) ) dt dt Fcp(TR - Tfl) - V∆h + Qfl (7)
Ki )
dH ) Fhf - Vhv - Lhl + Qfl dt
dH
Mfl cp
i
dt
(6)
where ni are the mole numbers of components A and B and (∂H/∂ni) are their partial molar enthalpies, which are equal to the pure component molar enthalpies for ideal mixtures to be considered in this paper. Thus, for constant pressure pfl and constant holdup Mfl, the energy balance finally reduces to
The saturation pressures pis are calculated by the ClausiusClapeyron relation
ln
(
pis ∆hi 1 1 ) pfl R Ti Tfl
)
(12)
Here, Ti is the boiling point temperature of component “i” at pressure pfl and pis is the saturation pressure at Tfl. The motivation behind using the Clausius-Clapeyron equation is that all parameters occurring in it have a direct physical meaning and thereby greatly help with the interpretation of the results. The set of eqs 1, 2, 7, 9, and 10 represents a system of differential algebraic equations (DAE), which has to be solved for the unknowns xe, ye, V, L, Tfl. In this model formulation, xe is a differential variable, which follows from the differential eq 1. For constant pressure, the temperature Tfl is an algebraic variable, which follows from the equilibrium condition (eq 10). The flow rates V and L are thus calculated from the total mass and energy balances (eqs 2 and 7). However, since the energy balance is a differential equation for xe and Tfl, the DAE system has a differential index of 2 (ref 8) and is therefore hard to solve numerically.9 There are two possibilities to overcome this index-2-problem. These are either a formal index reduction or a model reformulation by modifying the underlying assumptions. In the present case, a quasistatic energy balance, for example, will lead to a differential index of 1 (model 1, quasistatic energy balance). The quasistatic energy balance reads
0 ) Fcp(TR - Tfl) - V∆h + Qfl
(13)
For a formal index reduction (model 2, rigorous energy balance), the time derivative of the temperature in the energy balance (eq 7) is calculated by implicit differentiation of the equilibrium condition (eq 10). This gives rise to the following expression for the temperature derivative
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1199
dTfl dTfl dxe ) dt dxe dt )-
(14) RTfl2(KA - KB)
dxe ∆hAKAxe + ∆hBKB(1 - xe) dt
(15)
2.2. Stability. Let us first analyze the stability of the two models discussed above. Model 1: Quasistatic Energy Balance. The uniqueness and stability of the steady states of the constantly heated stand alone flash with quasistatic energy balance has been established in our previous work.1 Model 2: Rigorous Energy Balance. The stability analysis is now extended to the stand alone flash with dynamic energy balance. From the total material balance (eq 2), L can be calculated and inserted into eq 1. Then, the material balance for component A can be written as
Mfl
dxe ) F(z - xe) + V(xe - ye) dt
(16)
The energy balance (eq 7) can be rewritten as
(
Mfl cp
)
dTfl dxe + (hl,A - hl,B) ) dxe dt Fcp(TR - Tfl) - V∆h + Qfl (17)
From eqs 16 and 17, the unknowns xe and V are calculated. Tfl and ye are calculated separately from eqs 9 and 10 in dependence of xe. The stability of the stand alone flash with rigorous energy balance follows from the eigenvalues of the corresponding linearization of eqs 16 and 17. These are obtained from the generalized eigenvalue problem
det
[( ) ( )] f f aλ 0 - f11 f12 bλ 0 21 22
)0
b ) Mflcp
dTfl dxe
dye dxe
f21 ) -Fcp f22 ) -∆h
the partial derivatives of the right-hand sides with respect to xe and V. Through evaluation of eq 18, λ is obtained as
λ)
f11f22 - f12f21 af22 - bf12
which is always negative because of the stability of model 1 as discussed above. Further, it is easily shown that af22 is always negative. For a model with the dynamic energy balance (model 2), b is not equal to 0 and has to be taken into account in eq 25. However, since bf12 is also always positive, model 2 also always has a stable steady state. The positivity of bf12 follows from the following facts. For the binary nonazeotropic mixtures considered here, the product dTfl/dxe(xe - ye) is always positive. If component A is the heavy component, the temperature derivative dTfl/dxe and xe - ye are both positive (see Figure 2 in ref 1). If component A is the light boiling component, the temperature derivative dTfl/dxe and xe - ye are both negative. Further, it should be noted that the derivative dye/dxe is always positive for thermodynamically stable mixtures.1 2.3. Dynamic Transient Behavior. In this subsection, the dynamic transient behavior of the two models described above is considered. As a reference, a third model is introduced with variable pressure pfl inside the flash. Here, pfl is calculated from a pressure drop relation according to10
V2 V pfl - p0 ) 0.9∆pn 2 + 0.1∆pn D Dn n
(27)
3. Reactor-Separator System with Liquid Recycle
(20)
(22)
dye dTfl - V(∆h1 - ∆h2) dxe dxe
(26)
(24)
(19)
(21)
f12 ) xe - ye
f11f22 - f12f21 af22
(23)
(18)
of the left-hand sides and
f11 ) -F + V - V
λ)
where p0 is the constant external pressure, ∆pn is the standard pressure loss, and Dn is the standard throughput. Equation 8 is used as an energy balance. The introduction of the effect of pressure variation in fact makes the problem easier to solve numerically, since now the DAE system has differential index of 1. This can be seen readily since the pressure drop equation is used to determine V, the VLE is used to obtain pfl, and the temperature of the flash is obtained from the dynamic energy balance. The steady-state compositions predicted by the three models for the stand alone flash are the same. Moreover, the steady states are also found to be stable. Hence, any difference is going to be in the dynamic response of the various models to disturbances. To study this, the system is allowed to attain a steady state at a feed flow rate of F0 ) 25 mol/s (Figure 2). At t ) 15 s, a step input is applied to change the fresh feed flow rate to F0 ) 35 mol/s. It is seen that the quasistatic energy balance responds immediately to the disturbance while the dynamic energy balances take some time before attaining the new steady state. Similarly when the flow rate is changed back to the original value, it is found that the quasistatic energy balance is the one with the fastest response and the other dynamic balances are slower. In fact, the response of the dynamic energy balance with a constant pressure as well as with the pressure drop included are almost identical.
with the coefficients
a ) Mfl
For the quasistatic energy balance, b is equal to 0 and model 2 becomes equivalent to model 1. The eigenvalue λ reduces to
(25)
In this section, the focus is on a reactor-separator system with liquid recycle as illustrated in Figure 1b. The reactor is an isothermal continuous stirred tank reactor with a first-order chemical reaction A f B. The separator is a single stage constantly heated flash as in the previous section. F0 represents the fresh feed stream to the system. The holdup of the reactor is MR, and the mole fraction of component A in the tank is
1200
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007
Figure 2. Stand alone flash with different energy balances: (model 1) quasistatic energy balance; (model 2) rigorous energy balance with constant pressure; (model 3) rigorous energy balance with variable pressure. The parameters are given in Table 1. F ) 25.0 mol/s until t ) 15 s. From t ) 15-25 s, F ) 35.0 mol/s. At t ) 25 s, F ) 25.0 mol/s again. The solid line indicates the stable steady state. Table 1. Parameters Values Used for the Numerical Examples value param MR xf R TR ∆hA ∆hB TA TB T0 Mfl cp Qfl p0
Table 2. Physical Properties Used for Acetic Acid, Ethanol, and Water parameter
Fig 2
Fig 3,4
Fig 5,6
Fig 7,8
Fig 9,10, 12-14
unit
0.3 ) z 8.3144 280.0 100000.0 25080.0 400.0 300.0 273.15 20.0 209.836 15 × 105 101325.0
100.0 0.3 8.3144 385.0 25080.0 100000.0 400.0 300.0 273.15 20.0 209.836 15 × 105 101325.0
100.0 0.3 8.3144 310.0 100000.0 25080.0 400.0 300.0 273.15 20.0 209.836 15 × 105 101325.0
100.0 0.1 8.3144 310.0 100000.0 25080.0 400.0 300.0 273.15 20.0 209.836 15 × 105 101325.0
100.0 0.3 8.3144 310.0 25080.0 100000.0 400.0 300.0 273.15 20.0 209.836 15 × 105 101325.0
mol dimless J/(K mol) K J/mol J/mol K K K mol J/(K mol) W Pa
given by z. The molar flow rate leaving the tank is F. The molar liquid and vapor outflows leaving the flash are L and V. Component A is the heavy boiling component. The equilibrium concentration of component A in the liquid and the vapor of the flash unit are given by xe and ye, respectively. The liquid holdup of the flash is Mfl, and its pressure and temperature are given by pfl and Tfl. In a previous paper, it was shown that the instability and multiplicity of steady states of the coupled system in Figure 1 may arise because of the recycle if the reactor temperature TR is either lower than the flash temperature Tfl or the enthalpy of vaporization of the heavy boiling component is smaller than the enthalpy of vaporization of the light boiling component.1 The first condition, i.e., TR < Tfl, is always satisfied if the reactor is operated in the region of a single liquid phase under the same pressure as the separator, since the latter is in the two phase region. The second condition, i.e., ∆hA < ∆hB, is satisfied only in special cases. Mixtures of acetic acid and water as well as acetic acid and ethanol are typical examples to be treated subsequently (see also Table 2).
normal boiling point (K) enthalpy of vaporization (J/mol) at boiling point molar heat capacities (J/(mol K)) at 298.15 K
symbol
acetic acid
ethanol
water
Ti ∆hi
391.1 23700.0
351.4 38560.0
373.2 40650.0
cp,i
123.3
112.3
75.3
The analysis in ref 1 was based on a quasistatic energy balance. In the present paper, also a rigorous energy balance is applied and compared with the quasistatic energy balance. It will be shown that in contrast to the stand alone flash the dynamic behavior of the system may crucially depend on the specific model formulation. To alleviate the analysis and illustrate some basic phenomena, we focus on the limit of vanishing chemical reaction. Thus, the reactor serves as a pure mixing tank, where the liquid from the recycle is mixed with the fresh feed. The influence of finite reaction is discussed in the conclusion section. 3.1. Model Equations. Assuming isothermal operation at TR, the model equations of the mixing tank follow from the corresponding material balances. Material balance of component A
MR
dz ) F0xf + Lxe - Fz dt
(28)
Total material balance
0 ) F0 + L - F
(29)
From these equations, it is obvious that the stand alone mixing tank represents a linear system, which always admits a unique and stable steady state. The downstream flash unit is modeled using the same equations as described before. Once again, two different model formulations of the energy balance of the flash are considered as before to analyze the stability of the coupled system.
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1201
3.2. Stability/Dynamic Behavior. Model 1: Quasistatic Energy Balance. For fixed fresh feed flow rate F0 and constant holdups in the mixing tank and the flash drum, the vapor flow rate V has to be equal to the flow rate F0 at every instant of time. This follows from the overall material balance of the plant. F, the liquid flow rate leaving the mixing tank, can be obtained from the energy balance of the flash unit (eq 13) as
F)
F0∆h(ye(xe)) - Qfl
play an important role in a industrial applications like acetic acid production or esterification processes. Model 2: Rigorous Energy Balance. The system is now governed by three dynamic balances which consist of the material balance of the tank, the material balance of the flash, and the rigorous energy balance of the flash
MR
dz ) F0xf - Fz + (F - F0)xe dt
(39)
Mfl
dxe ) Fz - F0ye - (F - F0)xe dt
(40)
(30)
cp(TR - Tfl(xe))
The liquid recycle flow rate from the flash L is
L ) F - V ) F - F0
(31)
Using relation 5 and the equilibrium relations 9 and 10, the flow rates F and L can be expressed as functions of xe. In view of eqs 30 and 31 and V ) F0, the component material balances for the mixing tank (eq 28) and the flash (eq 1) can be written as a two-dimensional dynamic system
dz ) F0xf + (F(xe) - F0)xe - F(xe)z dt
(32)
dxe ) F(xe)z - F0ye(xe) - (F(xe) - F0)xe dt
(33)
MR Mfl
(
Stability follows from the Jacobian of eqs 32 and 33
-
A)
F MR
F Mfl
(
dF 1 L + (xe - z) MR dxe
(
)
dye 1 dF -L - F0 - (xe - z) Mfl dxe dxe
)
)
Here, the unknowns are z, xe, and F, where F is obtained from
Mfl cp
dTfl dxe ) Fcp(TR - Tfl) - F0∆h + Qfl dxe dt
the energy balance of the flash. Stability follows from the eigenvalues of the linearization of eqs 39, 40, and 41. These are obtained from the generalized eigenvalue problem
(
)
aλ - f11 -f12 -f13 -f bλ f det 21 22 -f23 ) 0 -f31 cλ - f32 -f33
a ) MR
(43)
b ) Mfl
(44)
(34) c ) Mfl cp
det(A) > 0
(35)
and trace(A) < 0
(36)
dTfl dxe
(45)
of the left-hand side and
The first condition ensures static stability. The determinant
F0F dye MRMfl dxe
(42)
with the coefficients
Necessary and sufficient conditions for the dynamic stability of a second-order system are (e.g., see ref 11 page 311)
det A )
(41)
f11 ) -F
(46)
f12 ) L
(47)
f13 ) xe - z
(48)
f21 ) F
(49) dye dxe
(37)
f22 ) -L - F0
is always positive. Therefore, the system is statically stable. Bifurcations leading to steady-state multiplicity, which have been found for the reactive case,1 can therefore not occur here. In addition to static stability, the trace of the Jacobian has to be negative for dynamic stability. This depends on the sign of the derivative dF/dxe. From eq 30, one obtains
f23 ) z - xe
(51)
f31 ) 0
(52)
(
)
dye dTfl 1 dF ) F0(∆hA - ∆hB) + Fcp dxe cp(TR - Tfl) dxe dxe
f32 ) -Fcp
f33 ) cp(TR - Tfl)
(38)
Hence, similar to the reactive case,1 dynamic stability can only be proven if TR > Tfl and ∆hA > ∆hB corresponding to dF/dxe > 0 in eq 38. On the reverse, dynamic instability in the form of self-sustained oscillations can occur if TR < Tfl or ∆hA < ∆hB. As discussed above, the first condition is always met if the reactor and separator are operated under the same pressure, whereas the second condition is only satisfied for specific mixtures, like acetic acid-water, for example, which however
dye dTfl - F0(∆hA - ∆hB) dxe dxe
(50)
(53) (54)
the partial derivatives of the right-hand sides with respect to xe, z, and F. This yields the following characteristic equation for the determination of the eigenvalues λ.
λ2x + λy + z ) 0
(55)
y ) af32f23 + c(f11f23 - f13f21) - f33(f11b + f22a)
(56)
with the coefficients
1202
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007
x ) a[bf33 - cf23] )
[
MR Mflcp(TR - Tfl) + Mflcp
]
dTfl (x - z) (57) dxe e
z ) f33f11f22 + f32f13f21 - f32f11f23 - f12f21f33
(58)
The solutions of these quadratic equations are
-y ( xy2 - 4xz λ1,2 ) 2x
(59)
Table 3. Location of Bifurcation Point for Different Energy Balancesa energy balance
fig
quasistatic rigorous (pfl ) const) quasistatic rigorous (pfl ) const) quasistatic rigorous (pfl ) const) quasistatic rigorous (pfl ) const)
3 4 5 6 7 8 9a 9b
F0 value of pole of eigenvalues (mol/s)
F0 value of Hopf point (mol/s)
19.1768
24.2580 24.2580 16.3957 16.3957 18.3024 18.3024 6.7992, 12.6093
20.4554 3.5095
a
For the special case of a quasistatic energy balance, c is equal to 0 and the characteristic equations of model 1 and 2 are equivalent. Thus, the system is stable if TR > Tfl and ∆hA > ∆hB. The eigenvalues in eq 55 have a pole for Tfl ) TR, which was therefore excluded so far. For the dynamic energy balance, c is not equal to 0. Additional poles of the eigenvalues in eq 55 with a change of sign will arise for
Mflcp(Tfl - TR) ) Mflcp
dTfl (x - z) dxe e
For the case of the rigorous energy balance with variable pressure, the same bifurcations have been found at the same values as in the case of the rigorous energy balance with constant pressure. ∆pn ) 100.0 Pa and Dn ) 25.0 mol/s have been used for the model with variable pressure.
(60)
Since the right-hand side of this equation is always positive for nonazeotropic mixtures as discussed above (see also Figure 2 in ref 1), such a pole will arise only for Tfl > TR, i.e., in cases where instability is possible. In the following, implications are illustrated with numerical examples. In particular, all combinations are considered where instability is possible, i.e., (i) Tfl < TR and ∆hA < ∆hB, (ii) Tfl > TR and ∆hA > ∆hB, and (iii) Tfl > TR and ∆hA < ∆hB. Especially, the role of poles of the eigenvalues for Tfl > TR will be studied. For each of these cases, the bifurcation behaviors predicted by a quasistatic energy balance and a dynamic energy balance at constant pressure are compared. It should be noted that a dynamic energy balance with variable pressure as introduced in the previous section for the stand alone flash gives almost identical solutions as a dynamic energy balance with constant pressure and is therefore not mentioned explicitly in the remainder. The parameter sets used for the simulations are given in Table 1 and the respective figure caption. The locations of the bifurcation points are summarized in Table 3. (i) Tfl < TR and ∆hA < ∆hB. First, instabilities arising for Tfl < TR and ∆hA < ∆hB are considered. Here, the instability arises due to the violation of the condition for the heats of vaporization. In this situation, the heat of vaporization of the heavy boiling component is less than that of the light boiling component. The dependence of the steady state and its stability on the fresh feed flow rate F0 is depicted in Figures 3 and 4. It is seen that under the conditions of a quasistatic energy balance the system admits a Hopf point where a subcritical bifurcation occurs at F0 ) 24.2580 mol/s (Figure 3). The steady state is stable for flow rates larger than this value and is unstable for flow rates lower than this value. The bifurcation behavior when using the dynamic energy balance is depicted in Figure 4. This also shows that the system exhibits a similar kind of behavior. The steady state is again unstable for F0 less than 24.2580 mol/s and stable for F0 more than 24.2580 mol/s. A subcritical Hopf bifurcation occurs at this point. For this case, it can be seen that the bifurcation characteristics are independent of the choice of the energy balance. Both models predict the same kind of behavior. Practical examples, where the enthalpy of vaporization of the heavy boiling component is smaller than the enthalpy of
Figure 3. Liquid recycling configuration: dependence of xe (black line), ye (red line), and z (green line) on F0. The parameters are given in Table 1: (solid line) stable steady state; (dashed line) unstable steady state; (square) Hopf bifurcation; (open circles) maxima of unstable periodic oscillations.
Figure 4. Liquid recycling configuration with rigorous energy balance (pfl ) const): dependence of xe (black line), ye (red line), and z (green line) on F0. The parameters are the same as in the case of Figure 3 (Table 1): (solid line) stable steady state; (dashed line) unstable steady state; (square) Hopf bifurcation; (open circles) maxima of unstable periodic oscillations. The branch of unstable oscillations (open circles) emerges at a subcritical Hopf bifurcation and crosses the feasibility boundary ye ) z at F0 ) 24.3807 mol/s.
vaporization of the light boiling component are mixtures of acetic acid and water or acetic acid and ethanol, respectively. In both cases, acetic acid is the heavy boiling component A.
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1203
Figure 5. Liquid recycling configuration: dependence of xe (black line), ye (red line), and z (green line) on F0. The parameters are given in Table 1: (solid line) stable steady state; (dashed line) unstable steady state; (square) Hopf bifurcation; (asterisks) maxima of stable periodic oscillations; (circles) maxima of unstable periodic oscillations. The branch of stable oscillations (asterisks) crosses the feasibility boundary ye ) z at F0 ) 16.345 mol/s.
The relevant physical properties for these components were taken from ref 12 and are summarized in Table 2. The heat capacity of the mixture was assumed to be constant. The arithmetic mean of the individual heat capacities of the components was taken for the calculations. For the acetic acidwater system, a total heat capacity cp of 99.3 J/(K mol) was taken, and for the acetic acid-ethanol system, a heat capacity of 117.8 J/(K mol) was taken. The mole fraction xf of compound A (acetic acid) in the feed was set to 0.5. All other parameters were the same as for Figures 3 and 4 (see Table 1). In both cases, qualitatively similar bifurcation diagrams as discussed above are obtained. Hopf bifurcation points were found at F0 ) 47.7403 mol/s for the acetic acid-water system and at F0 ) 51.2435 mol/s for the acetic acid-ethanol system for all three models. (ii) Tfl > TR and ∆hA > ∆hB. Next, the case is discussed when only the condition on the temperature is violated. Here, Tfl > TR and ∆hA > ∆hB is taken. The corresponding bifurcation diagrams are shown in Figures 5 and 6. For the quasistatic energy balance, the bifurcation diagram exhibits a subcritical Hopf bifurcation and a branch of unstable oscillations emerges (Figure 5). This branch of limit cycles exhibits a cyclic fold. In this oscillatory state, the branch of oscillations crosses the boundary z ) ye at F0 ) 16.345 mol/s where it becomes infeasible. For the same operating conditions in the case of the dynamic energy balance, the system exhibits a Hopf bifurcation at F0 )16.3957 mol/s as well as a pole of the eigenvalues at 19.1768 (Figure 6). These diagrams represent the behavior of the system when the feed composition xf equals 0.3. Now, the mole fraction of A in the feed xf is reduced to 0.1. In this case, the bifurcation diagram exhibits a supercritical Hopf point at which a branch of stable periodic solutions emerges (Figure 7). This branch seems to end at a homoclinic bifurcation. This is supported by the fact that the oscillation period increases strongly while following the branch of oscillations as the bifurcation parameter is increased beyond the bifurcation point. In contrast to that, a subcritical Hopf bifurcation at F0 ) 18.3024 mol/s and a pole of the eigenvalues at F0 ) 20.4554 are obtained for the case of a dynamic energy balance as shown in Figure 8. (iii) Tfl > TR and ∆hA < ∆hB. Next, the case is considered when both conditions are violated, i.e., Tfl > TR and ∆hA
TR holds. In that case, a Hopf bifurcation is found at F0 ) 52.691 mol/s for the case of the quasistatic energy balance (model 1). For the case of the rigorous energy balance (model 2), both a Hopf bifurcation at F0 ) 52.691 mol/s and a pole of eigenvalues with a sign change at F0 ) 53.04 mol/s emerge. A stable steady state is found at F0 values below the Hopf bifurcation, whereas an unstable steady state is found at higher F0 values. At F0 ) 53.04 mol/s, no stability change occurs.
So far, transportation delays between the different units in Figure 1b were neglected. In previous work,4,5 the influence of finite transportation delay for a system similar to Figure 1b was studied. In addition, a chemical reaction in the mixing tank was taken into account.4,5 In particular, an isothermal first-order reaction, an isothermal cubic autocatalytic reaction,4 and an exothermic reaction5 were considered. Isothermal and isobaric operation of the flash was assumed and a delay was considered in the forward branch between the reactor and the separator. In contrast to this, the focus here is on the nonreactive case, with a constantly heated flash, and a delay in the backward branch, i.e., between the flash and the reactor, is considered. A finite delay in the coupled system arises since the transport velocity of the streams moving from one unit to another is finite. This delay will not affect the steady-state solution values and the steady-state branches. Although the stationary values of the state variables remain the same, the delay may change the stability of the steady state. Dynamic instability can arise in a system when the delay is above a certain critical value in a coupled plant. This will also depend on the flow control strategies of the plant (see also ref 1). As shown in the previous section, the nonreactive liquid recycle system illustrated in Figure 1b can show self-sustained dynamic behavior in the absence of delays. Now, the behavior of the system in the presence of a delay in the recycle flow is analyzed. This implies that changes in the separation unit do not immediately enter the mixing tank but appear only after a delay of some time. The delay time is directly proportional to the length of the tube in Figure 11. The corresponding model equations are
dz ) F0xf + (F - F0)xe(t - τ) - Fz dt
(61)
dxe(t) ) Fz - F0ye(xe(t)) - (F - F0)xe(t) dt
(62)
MR Mfl
where τ is the delay time. It should be noted that for an incompressible liquid only the concentrations in the liquid recycle are delayed. The flow rate at the outlet of the recycle
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1205
Figure 11. Schematic representation of the liquid recycle system with finite transportation delay.
tube will respond instantaneously to changes of the flow rate at the inlet of the tube. Hence, for a finite delay, the flow rate F in eqs 61 and 62 follows from eq 38 and is a function of the instantaneous value xe(t). The characteristic equation of a delay system is13
(∂x∂f + ∂y∂f e
det
-λτ
)
- λI ) 0
(63)
Here, ∂f/∂x represents the Jacobian with respect to the undelayed variables and ∂f/∂y represents the Jacobian with respect to the delayed variables. I is the identity matrix. With respect to eqs 61 and 62, it is found that
(
F -λ MR det F Mfl -
(
)
)
1 dF Le-λτ + (xe - z) MR dxe )0 dye 1 dF -L - F0 - (xe - z) -λ Mfl dxe dxe (64)
(
)
which is almost identical to eq 34 except for the factor e-λτ in the first row. This represents a transcendental equation for the eigenvalues λ. The solution set of this characteristic equation is a whole spectrum of eigenvalues. Stability conditions (eqs 35 and 36) only apply to a quadratic characteristic equation and are therefore no longer valid. For the delay system to be stable, the whole spectrum of eigenvalues has to lie in the left half plane. Since an analytical approach seems not possible anymore, in the remainder, a numerical approach is used. For the numerical solution of the delay differential equations (eqs 61 and 62), different approaches are possible. This could be either a direct approach using specialized software for delay differential equation systems like the dde23 package for MATLAB,14 for example, or an indirect approach using standard software for ordinary differential equations. In this paper, the latter approach is used. For this purpose, the transport in the tube in Figure 11 is modeled with a partial differential equation according to
∂x ∂x +V )0 ∂t ∂z
(65)
5. Discussion and Concluding Remarks
(66)
In this paper, the role of the energy balance formulation on the dynamic behavior of evaporative separators was studied. The steady state behavior predicted by the different models was all the same. It was found that the specific formulation of the
with boundary values
x(z ) 0, t) ) xe(t), x(z ) l, t) ) xe(t - τ)
In particular, a transport velocity V of 1 m/s is chosen. Then, the delay time is directly equal to the length of the tube. For the numerical solution, the partial differential equation eq 65 is converted to a set of ordinary differential equations using a method of lines. A finite volume method with total variation diminishing (TVD) was applied15 for obtaining the numerical solution. The number of grid points was set to nz ) 400. Convergence of the results was checked by comparing these results with those obtained using 800 grid points, which yields almost the same results. The numerical simulation of the resulting set of ordinary differential equations together with eqs 61 and 62 was accomplished with DIVA.16 Now, the behavior of the liquid recycle system with finite delay is discussed. First, we assume a quasistatic energy balance. To investigate the influence of finite delay, the parameter set of Figure 9 is used. A value of F0 ) 8.0 mol/s is chosen where periodic oscillations occur when the system has no delay (see Figure 9). The delay time τ is increased successively. The resulting dynamic behavior is calculated by time integration of the discretized system as described above. In Figure 12, the reduced period (oscillation period with delay P divided by period without delay P0) is depicted versus the delay time τ. The reduced periods increase smoothly and decrease sharply in a sawtoothlike manner. By stepwise increasing and decreasing the parameter τ, regions were found where the different branches in Figure 12 overlap. Within these regions, two or even more stable limit cycles coexist for the same set of parameters. This phenomenon is called bi- or multirhythmicity, respectively, and is the analogy to bi- or multistability of steady states. Figure 13a shows the complete bifurcation diagram. A branch of stationary solutions occurs at Tfl ) 345.08 K. These stationary solutions are only stable between the Hopf points at τ ) 0.22 and 0.48 s. Otherwise, one or more stable periodic solutions exist in addition to the unstable steady state. For τ ) 0.220.48 s, the oscillations vanish and give way to a stable stationary solution. The Hopf bifurcation points at τ ) 0.48 and 2.12 s as well as the Hopf bifurcation points at τ ) 1.96 and 3.68 s are connected by a branch of periodic oscillations. At the Hopf point at τ ) 3.66 s, another branch of oscillations emerges. The birhythmic and multirhythmic regions are indicated by the overlapping branches which indicate the maxima of temperature of the stable periodic oscillations. These branches are connected by branches of unstable periodic oscillations (not shown). Figure 14 elucidates the birhythmicity in the system. It shows switching between two different stable periodic solutions for a value of τ ) 2.05 s for the case of the quasistatic energy balance. At t ) 100 s, the initial conditions of the variables are changed so that a transition to the other limit cycle takes place. Switching back to the original initial values would again cause the system to switch to the original oscillatory state. Thus, it can be seen that delay can induce birhythmicity in the system. In the next step, the rigorous energy balance (pfl ) const) is applied. Figure 13b shows the bifurcation diagram for the case of a rigorous energy balance using the same parameters as in Figure 13a. Without delay, an unstable steady state exists. By increasing the delay, the unstable steady state remains unstable. The system is monotonically unstable, and no further attracting solutions have been found by dynamic integration.
1206
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007
Figure 12. Liquid recycling configuration with delay: dependence of the reduced oscillation period (P/P0) on the delay time τ. F0 ) 8.0 mol/s. The oscillation period without delay is P0 ) 2.11 s. All other parameters are as in Table 1. In the regions where different branches are overlapping, birhythmicity or multirhythmicity occurs.
Figure 14. Liquid recycling configuration with delay and a quasistatic energy balance (model 1): switching between two stable oscillatory states in the birhythmic region at t ) 100 s by using a different set of initial conditions for the variables. F0 ) 8.0 mol/s, and τ ) 2.05 s. All other parameters are as given in Table 1.
Figure 13. (a) Liquid recycling configuration with delay and a quasistatic energy balance (model 1). (b) Liquid recycling configuration with delay and a rigorous energy balance (model 2): dependence of Tfl versus the delay time τ. F0 ) 8.0 mol/s. All other parameters are as given in Table 1: (solid line) stable steady state; (dashed line) unstable steady state; (asterisks) maxima of stable periodic oscillations.
energy balance of such a system plays a minor role if a stable stand alone flash is considered. In contrast to this, it can play an important role if potentially unstable multiunit plants with recycles are considered. Some deeper insight was obtained from a rigorous analysis of a suitably simplified model system which was a mixing tank with a single stage flash and recycle. From the analysis, we were able to derive a condition, under which different results have to be expected. The different behavior results from a pole of the eigenvalues with a change of sign. In tendency, a system with a rigorous dynamic energy balance turned out to be most of all monotonically unstable, whereas a simplified model with a quasistatic energy balance more often admits “bounded” instability in the form of oscillatory solutions.
From the present study, two important conclusions can be drawn. First, the various model formulations discussed in this paper directly carry over to multistage distillation processes, which are of major importance in practice. Similar phenomena, reported here for simplicity for a single stage flash process, can also be found if a multistage distillation process is used instead. In particular, it was found that a model with a quasistatic energy balance can show complex periodic and even chaotic behavior, whereas a model with a rigorous dynamic energy balance is monotonically unstable. Similar patterns of behavior were observed if in addition an isothermal first-order reaction in the mixing tank is taken into account. The details of this will be given elsewhere.17 Second, models which give realistic predictions for the simulation of stand alone process units may fail for the simulation of multiunit processes with recycles. Hence, for dynamic flowsheet simulation, further model validation is required. Acknowledgment The authors would like to thank the Volkswagen-Stiftung for financial support under grant I/77311. A.K. acknowledges the financial support of the Alexander-von-Humboldt Foundation.
Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1207
Notation
Literature Cited
A ) Jacobian matrix cp ) molar heat capacity (J/K/mol) Dn ) standard throughput (mol/s) F ) feed to the flash (mol/s) F0 ) feed to the mixing tank (mol/s) H ) total enthalpy (J) h ) molar enthalpy (J/mol) ∆h ) enthalpy of vaporization (J/mol) I ) identity matrix Ki ) equilibrium constant (-) L ) flow rate of the liquid from the flash (mol/s) l ) length (m) λ ) eigenvalue (1/s) M ) molar holdup (mol) n ) number of moles (mol) nz ) number of grid points (-) pfl ) pressure inside the flash (Pa) pis ) vapor (saturation) pressure of pure component “i” (Pa) p0 ) external pressure (Pa) ∆pn ) standard pressure loss (Pa) Q ) heating rate (W) R ) gas constant (J/(K mol)) T ) temperature (K) Ti ) boiling point temperature of component “i” (K) T0 ) reference temperature (K) t ) time (s) τ ) delay time (s) V ) flow rate of the vapor from the flash (mol/s) V ) flow velocity (m/s) xe ) mole fraction of reactant A in the liquid from the flash (-) xf ) mole fraction of reactant A in the feed (-) ye ) mole fraction of reactant A in the vapor from the flash (-) z ) mole fraction of reactant A in the liquid from the mixing tank (-) z ) spacial coordinate, eqs 65 and 66 (m)
(1) Zeyer, K.-P.; Pushpavanam, S.; Kienle, A. Nonlinear behavior of reactor-separator networks: influence of separator control structure. Ind. Eng. Chem. Res. 2003, 42, 3294. (2) Pushpavanam, S.; Kienle, A. Nonlinear behavior of an ideal reactor separator network with mass recycle. Chem. Eng. Sci. 2001, 56, 2837. (3) Sagale, A. A.; Pushpavanam, S. A comparison of control strategies for a nonlinear reactor-separator network sustaining an autocatalytic isothermal reaction. Ind. Eng. Chem. Res. 2002, 41, 2005. (4) Balasubramanian, P.; Kosuri, M. R.; Pushpavanam, S.; Kienle, A. Effect of delay on the stability of a coupled reactor-separator system. Ind. Eng. Chem. Res. 2003, 42, 3758. (5) Balasubramanian, P.; Pushpavanam, S.; Kienle, A.; Balaraman, K. S. Effect of delay on the stability of a coupled reactor-flash system sustaining an elementary non-isothermal reaction. Ind. Eng. Chem. Res. 2005, 44, 3619. (6) Krishna, M. V.; Pushpavanam, S.; Kienle, A.; Vetukuri, S. R. R. Nonlinear behavior of reactor-separator systems with azeotropic mixtures. Ind. Eng. Chem. Res. 2006, 45, 212. (7) Vetukuri, S. R. R.; Pushpavanam, S.; Zeyer, K.-P.; Kienle, A. Nonlinear behavior of coupled reactor-separator systems with azeotropic vapor-liquid equilibriums (VLEs): comparison of different control strategies. Ind. Eng. Chem. Res. 2006, 45, 1019. (8) Unger, J.; Kro¨ner, A.; Marquardt, W. Structural analysis of differential-algebraic equation systems: theory and applications. Comput. Chem. Eng. 1995, 19, 867. (9) Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solution of Initial Value Problems in Differential-Algebraic Equations; North Holland & Elsevier Science Publishing Company: New York, 1989. (10) Holl, P. Entwicklung und Einsatz eines dynamischen Simulators fu¨r Verfahrenstechnische Prozesse, Reihe 3; Verfahrenstechnik, no. 359, VDI-Verlag: Du¨sseldorf, 1994. (11) Bequette, B. W. Process Dynamics - Modeling, Analysis and Simulation; Prentice Hall PTR: Englewood Cliffs, NJ, 1998. (12) Handbook of Chemistry and Physics, 82nd ed.; CRC Press: Boca Raton, FL, 2001. (13) Epstein, I. R.; Pojman. J. A. An Introduction to Nonlinear Chemical Systems - Oscillations, WaVes, Patterns and Chaos; Oxford University Press: New York, 1998. (14) Shampine, L. F.; Thompson, S. Solving DDEs in MATLAB. Appl. Numer. Math. 2001, 37, 441. (15) Finlayson, B. A. Numerical Methods for Problems with MoVing Fronts; Ravenna Park Pub.: Seattle, WA, 1992. (16) Mangold, M.; Kienle, A.; Mohl, K. D.; Gilles, E. D. Nonlinear computation in DIVA- methods and applications. Chem. Eng. Sci. 2000, 55, 441. (17) Zeyer, K.-P.; Kulkarni, A. A.; Kienle, A.; Kumar, M. V.; Pushpavanam, S. Nonlinear behavior of reactor-separator and reactor-distillation networks:influence of the energy balance formulation. Submitted to the 17th European Symposium on Computer Aided Process Engineering (ESCAPE 17), Bucharest, Romania, May 27-30, 2007.
Subscripts f ) feed fl ) flash i ) component l ) liquid phase R ) mixing tank v ) vapor phase
ReceiVed for reView June 14, 2006 ReVised manuscript receiVed September 26, 2006 Accepted November 13, 2006 IE0607595
