Numerical Bifurcation Analysis of Delay in a Coupled Reactor

3920

Ind. Eng. Chem. Res. 2008, 47, 3920–3929

Numerical Bifurcation Analysis of Delay in a Coupled Reactor Separator System P. Balasubramanian* Department of Chemical Engineering, Indian Institute of Technology Guwahati, Guwahati 781 039, India

The effect of delay on the stability of a coupled reactor separator system is investigated in this work, using the DDE-BIFTOOL software package. The reactor considered is a continuously stirred tank reactor (CSTR), and the separator is modeled as a flash unit operating isothermally and isobarically. A first-order isothermal irreversible reaction A f B is considered to investigate the effect of delay arising as a result of transportation lag from the CSTR to the flash. The dynamic behavior of the coupled CSTR-flash system is governed by the delay differential equations. In this work, the dynamic instability induced by delay that is crucially dependent on the flow control strategy of the system is investigated using DDE-BIFTOOL. It is shown that delay can destabilize the steady state of the coupled system when the reactor effluent flow rate F is flow-controlled, and the molar holdup MR in the reactor is allowed to vary. When the fresh feed flow rate F0 is flow-controlled, the system exhibits delay-independent stability. The critical delay and frequency of oscillation obtained using DDE-BIFTOOL show good agreement with the same obtained using the analytical expressions proposed in the literature. 1. Introduction In a typical chemical process plant, a downstream separator is coupled with an upstream reactor via a reactant-rich recycle stream. Changes in the downstream separator can affect the upstream reactor through this recycle stream. The steady-state and dynamic behavior of a coupled reactor-separator system was investigated by several researchers. Luyben1 has analyzed the behavior of an isothermal continuously stirred tank reactor (CSTR) sustaining a first-order irreversible reaction, coupled with an isothermal-isobaric flash. It was found that the recycle flow rate was very sensitive to changes when the fresh feed flow rate (F0) is flow-controlled. Pushpavanam and Kienle2 have analyzed the behavior of a first-order exothermic reaction in a CSTR coupled with an isothermal-isobaric flash. They found that the coupled system could exhibit a maximum of two steady states when F0 is flow-controlled. Kiss et al.3 have investigated steady-state multiplicity in the CSTR-separator-recycle polymerization system sustaining a first-order isothermal reaction. It was observed that the feasible operation of the system was possible when the Damko¨hler number (Da) is greater than unity. Wu and Yu4 noted that the snowball effect does not disappear, but sometimes does manifest as large changes in the reactor holdup. A balance control structure was proposed in which the disturbance rejection effort was equally distributed between the reactor and the separator. The dynamic analysis and control of a reactor-distillation column with recycle was investigated by Kumar and Daoutidis.5 They proposed model reduction methodology for slow dynamics induced by large recycle streams. Bildea and Dimian6 analyzed the nonlinear behavior of several recycle systems involving firstand second-order reactions. They demonstrated that the reactor inlet was the appropriate location for fixing the flow rates in a coupled reactor-distillation-recycle system. This location ensures stable behavior and avoids the snowball effect. Seki and Naka7 proposed a hierarchical control structure for the reactor-separator-recycle system that was comprised of a lower-level regulatory control and a higher-level coordination control structure. The dynamic interaction between the reactor * To whom correspondence should be addressed. Tel.: +91-3612582263. Fax: +91-361-2690762. E-mail: [email protected].

and the separator can be minimized by fixing either the recycle flow rate or the reactor effluent flow rate in a regulatory control structure. Recently, a novel methodology was proposed for the design and control of recycle systems using nonlinear analysis.8 A clear distinction was made between self-regulation and feedback control regulation structures. It was shown that the best strategy for the plantwide control is fixing the reactor inlets on flow control. In the earlier studies, the effect of delay on the steady-state and dynamic behavior of a coupled reactor-separator system was neglected. Denn and Lavie9 have considered the effect of delay in the recycle loop of a coupled reactor-splitter system. They used linear transfer function models and focused on how the coupling affected the characteristic time and gain of the system. Balasubramanian et al.10 have studied the effect of delay on the stability of an isothermal CSTR sustaining an irreversible firstorder reaction A f B, coupled with an isothermal, isobaric flash. Here, composition delay was considered between the CSTR to the flash, and it was determined that the delay does not destabilize the steady state of the coupled system when the fresh feed flow rate F0 is flow-controlled. Also, it was observed that the delay can induce dynamic instability for other two control strategies: (i) when the reactor effluent flow rate (F) is flowcontrolled, and (ii) when the molar holdup in the reactor is allowed to vary. In the present work, DDE-BIFTOOL,11 which is a numerical toolbox for continuation and stability analysis of delay differential equations, is used to investigate the effect of delay which arises as a result of transportation lag from the reactor to the separator in a coupled reactor-separator system, as shown in Figure 1. As a result of delay, the composition entering the separator at any given time instant is that which was in the reactor at an earlier time instant. Hence, the recycle flow rate at a given time instant is a function of the reactor composition at an earlier time instant. The effect of delay arising due to the finite velocity of the recycle stream is neglected, as a result of fixed liquid composition leaving from the flash to the reactor instantaneously. The three control strategies that are considered to operate a coupled reactor-separator system are described as follows:

10.1021/ie071010c CCC: $40.75  2008 American Chemical Society Published on Web 05/02/2008

Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008 3921

The reactor considered is an isothermal CSTR, and the separator is a flash unit operated under isothermal and isobaric conditions. A first- order irreversible reaction A f B is considered to investigate the stability of the delay system. For isothermal operation of the CSTR, the coupled system does not show any dynamic instability in the absence of delay for the three control strategies depicted in Figure 1. The objectives of this work are (i) to investigate the stability of the delay system using DDE-BIFTOOL, and (ii) to compare the critical delay and frequency of oscillation obtained using DDE-BIFTOOL with the same obtained from the analytical expressions proposed in the literature.10 2. Theoretical Models In this section, the model equation that governs the behavior of the coupled reactor separator system is discussed. 2.1. Reactor Balance. The mole balance equation for the molar holdup MR, and the mole fraction of reactant A (z) in the CSTR are given as follows: dMR ) F0 + L - F dt′

Overall mole balance : Reactant A mole balance :

Figure 1. Typical flow control strategies of a coupled reactor-separator system: (a) F0 is flow-controlled, (b) F is flow-controlled, and (c) the variable reactor molar holdup (MR) case.

(1)

d(MRz) ) F0xaf + Lxe dt′ Fz(t ′ ) - kz(t ′ )MR (2)

2.2. Separator Balance. The flash is operated under isothermal and isobaric conditions, and, for a binary mixture, this uniquely fixes the composition of vapor and liquid. The molar holdup of the vapor and liquid phase in the flash is neglected and the two streams leaving the flash are assumed to be in equilibrium. Therefore, the overall mole balance and the mole fraction of reactant A balance of the flash are detailed below. F)V+L

(3)

Fz(t ′ -τ ′ ) ) Vye + Lxe

(4)

The transportation lag or delay from the reactor to the separator is considered in eq 4, and the effect of delay arising due to the finite velocity of the recycle stream is neglected in this analysis. When the reactor molar holdup MR is a constant, then the overall mole balance of the system becomes F0 ) V

Figure 2. Dependency of z on Da obtained using DDE-BIFTOOL for the fixed F, MR case (ye ) 0.2, xe ) 0.4, xaf ) 0.9, and Da ) 0.5).

(i) Fixed F0 and MR: In this case, the fresh feed flow rate (F0) is flow-controlled and the reactor effluent flow rate (F) is used as a manipulated variable to control the reactor molar holdup (MR) (see Figure 1a). (ii) Fixed F and MR: In this mode of operation, F is flowcontrolled, using a flow controller, and MR is controlled using F0 as the manipulated variable (see Figure 1b). (iii) Fixed F and F0: F0 and F are both flow-controlled (see Figure 1c). This fixes the recycle flow rate (L) at steady state, and the molar holdup MR of the reactor is determined by the interaction between the reactor and the separator.

(5)

In this work, it is assumed that ye < z < xe. This ensures that the flash can separate its feed stream into two streams of different compositions. The dimensionless delay differential equation that governs the dynamic behavior of the coupled CSTR-flash system operated under three different control strategies is described as follows. 2.3. Fixed F0, MR Case. The time evolution of the mole fraction of reactant A in the reactor is governed by the singlestate variable dimensionless delay differential equation when F0 is flow-controlled. It is represented as

(

) [(

z(t - τ) - ye dz(t) ) xaf + x dt xe - z(t - τ) e

) ]

xe - ye + Da z(t) xe - z(t - τ) (6)

where Da ) MRk/F0 and t ) t′F0/MR. The aforementioned equation is derived by eliminating the state variables F, L, and V using eqs 2–5 by substituting for them in terms of xe, ye, z, and F0.

3922 Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008

Figure 3. Roots of the characteristic equation (left) and real parts of the roots of the characteristic equation versus Da (right) plots for the fixed F, MR case obtained using DDE-BIFTOOL (ye ) 0.2, xe ) 0.4, xaf ) 0.9, and Da ) 0.5): (a and b) τ ) 1 and (c and d) τ ) 2.

2.4. Fixed F, MR Case. The governing single-state variable dimensionless delay differential equation of the mole fraction of reactant A in the reactor is

(

) (

)

xe - z(t - τ) z(t - τ) - ye dz(t) ) xaf + xe - (1 + Da)z(t) dt xe - ye xe - ye (7) where Da ) MRk/F and t ) t′F/MR. Equation 7 arises when we eliminate the state variables F0, L, and V using eqs 2–5 by substituting for them in terms of xe, ye, z, and F. 2.5. Fixed F0, F Case with Variable Molar Holdup. The evolution of the molar holdup (MR) and the mole fraction of reactant A (z) in the reactor are governed by the following multivariable dimensionless delay differential equations: dM/R dt

( ) ( )(

(

)1+r

z(t - τ) - xe xe - ye

)

) [( ) ) ] ( )(

xaf r z(t - τ) - ye 1 dz(t) ) + xe + / dt xe - ye MR M/R M/R r z(t - τ) - ye + Da z(t) xe - ye M/ R

(8)

The dimensionless variables are defined as Da ) MChk/F0, r ) F/F0, t ) t′F0/MCh, and MR* ) MR/MCh. The aforementioned equations are derived by eliminating the state variables L and V from eqs 1–4 by substituting for them in terms of xe, ye, z, F0, and F. 3. Stability Analysis of Delay Differential Equations Normally, delay differential equations have an infinite number of dimensions and the stability analysis of single-variable delay differential equations and multivariable delay differential equation is discussed in this section. 3.1. Single-Variable System. Consider the single-statevariable delay differential equation with a single constant delay of the form dz(t) ) f(z(t), z(t - τ), ψ) (10) dt subject to the condition z(t) ) z0 for -τ < t < 0. The characteristic roots of the delay differential equation are governed by the transcendental equation,12,13 ∂f ∂f + exp(-λτ) (11) ∂z(t) ∂z(t - τ) In the absence of delay, the steady state of the system is stable and the conditions under which delay can destabilize the steady λ)

(9)

Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008 3923

Figure 4. Real parts of the roots of the characteristic equation versus Da plots for the fixed F, MR case obtained using DDE-BIFTOOL (ye ) 0.2, xe ) 0.4, xaf ) 0.9, and Da ) 0.5): (a) τ ) 3, (b) τ ) 4, (c) τ ) 5, (d) τ ) 6, (e) τ ) 7, and (f) τ ) 8. I ) first branch, II ) second branch, and III ) third branch.

state of the system are determined. Bifurcations occur whenever the roots move through the imaginary axis as one or more parameters are changed. The system cannot be destabilized by a real λ crossing the imaginary axis (fold bifurcation), and, hence, delay can induce a Hopf bifurcation (i.e., the steady state becomes unstable when λ is purely imaginary). The critical delay

is obtained by equating the real and imaginary parts of the characteristic equation (eq 11) to zero. It is represented as ∂f ∂f + cos(ωτ) ∂z(t) ∂z(t - τ) ∂f sin(ωτ) 0)ω+ ∂z(t - τ)

0)

(12) (13)

3924 Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008

Figure 5. Dependency of dimensionless critical delay on Da (left) and frequency of oscillation versus dimensionless critical delay (right) for the fixed F, MR case (ye ) 0.2, xe ) 0.4, xaf ) 0.9, and Da ) 0.5): (a and d) τ ) 2-4, (b and e) τ ) 5-7, and (c and f) τ ) 8. Solid line represents data obtained using DDE-BIFTOOL, and symbol denote data obtained using eqs 21 and 22. I ) first branch, II ) second branch, and III ) third branch.

The frequency of oscillation ω and the critical delay τ determined explicitly from eqs 12 and 13 are ω)

∂f )
( ∂z(t∂f- τ) ) - (∂z(t)

(14)

) ⁄ ( ∂z(t∂f- τ) )) ∂f )
( ∂z(t∂f- τ) ) - (∂z(t)

(15)

2

((

cos -1 τ)

2

∂f ∂z(t)

2

2

Therefore, delay can destabilize the system when the argument of the arc cosine function lies in the interval from -1 to

+1 and ω > 0. The aforementioned stability analysis is applicable for the single-variable delay differential equation when F0 is flow-controlled, and F has flow-controlled strategies. 3.2. Multivariable System. Consider the system of delay differential equations with constant delays dzi(t) ) f(zi(t), zi(t - τj), ψ) (16a) dt where i varies from 1 to n and j varies from 1 to m. The roots of a characteristic equation that determines the stability of the steady state of the multivariable delay differential equations11–14 with delays is given by

Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008 3925

Figure 7. Dependency of z on r obtained using DDE-BIFTOOL for the variable reactor molar holdup case (ye ) 0.2, xe ) 0.4, xaf ) 0.9, Da ) 0.5, and r ) 1.5).

on the composition of the vapor and liquid streams leaving the flash, the delay-dependent effluent composition (z) from the CSTR to the flash, and the fresh feed flow rate F0 of the coupled system. It is represented as

(

L ) F0

z(t ′ -τ ′ ) - ye xe - z(t ′ -τ ′ )

)

(17)

The frequency of oscillation and critical delay are determined from the characteristic roots of the transcendental equation.10 It is given by ω)


( ) [( ) ] xe - ye xe - z

2

-

xe - ye + Da xe - z

2

(18)

and Figure 6. Mole fraction of reactant A (z) versus time (t) plots, showing the evolution to sustained oscillations for the fixed F, MR case (ye ) 0.2, xe ) 0.4, xaf ) 0.9, and Da ) 0.5): (a) τ ) 2.5 and (b) τ ) 6.

[

Det λI - J -

m

∑J j)1

τj

]

exp(-λτj) ) 0

τ) (16b)

The roots of the characteristic equation of the two statevariable delay differential equations with a single constant delay is determined by substituting m ) 1 in eq 16b. This step is essential to determine the roots of the characteristic equation of the delay differential equations that arise for the variable molar holdup case. 3.3. DDE-BIFTOOL. The DDE-BIFTOOL software package is a collection of Matlab routines for the numerical bifurcation analysis of systems of delay differential equations with several constant and state-dependent delays.11,14–16 The package implements the continuation and stability analysis of steady-state solutions, as well as fold and Hopf bifurcations. It further allows one to switch from the latter to an emanating branch of periodic solutions. 4. Results and Discussion In the following discussion, the stability analysis of the delay differential equation that governs the behavior of the coupled CSTR-flash system for the three control strategies is discussed. 4.1. Fixed F0, MR Case. In this mode of operation, the recycle flow rate (L) from the flash to the CSTR is dependent

{[( ) ] ⁄ ( )}
( ) [( ) ]

cos -1

xe - ye + Da xe - z

xe - ye xe - z

2

-

xe - ye xe - z

xe - ye + Da xe - z

(19)

2

In this control strategy, the frequency of oscillation is always a complex number for all positive values of Da, and the argument of the arc cosine is greater than unity. Hence, the system does not have a critical delay. Furthermore, at steady state, from eq 17, the recycle flow rate L lies in the interval from 0 to ∞ when the mole fraction of reactant A in the reactor approaches the bounds imposed by ye and xe from the flash, i.e., 0 < L/F0 < ∞ when ye < z < xe. Therefore, the coupled system can be operated within the range of complete conversion (z ) 0) and the zero conversion (z ) 1). Either increase or decrease in the fresh feed flow rate F0 entering into the reactor significantly affects the conversion of reactant A in the reactor. Also, L is highly sensitive to changes in the fresh feed flow rate (F0). Therefore, the effect of delay from the CSTR to the flash is insignificant, and, hence, this control strategy always exhibits delay-independent stability. 4.2. Fixed F, MR Case. In this case, the recycle flow rate L from the flash to the CSTR is given by

(

L)F

z(t ′ -τ ′ ) - ye xe - ye

)

(20)

3926 Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008

Figure 8. Real parts of the roots of the characteristic equation versus r plots for the variable reactor molar holdup case obtained using DDE-BIFTOOL (ye ) 0.2, xe ) 0.4, xaf ) 0.9, Da ) 0.5, and r ) 1.5): (a) τ ) 1, (b) τ ) 2, (c) τ ) 5, and (d) τ ) 12. I ) first branch and II ) second branch.

Here, the ratio of recycle flow rate L to the reactor effluent flow rate F lies in the interval from 0 to 1, when the mole fraction of reactant A in the reactor approaches the bounds imposed by ye and xe, i.e., 0 < L/F < 1, when ye < z < xe. On one hand, complete conversion of reactant A occurs in the reactor when L ) 0; on the other hand, zero conversion of reactant A occurs in the reactor when L ) F. Thus, feasible operation of this system is possible when L< F. The frequency of oscillation and the critical delay are found10 as ω) and


(

xe - xaf xe - ye

[(

cos -1 τ)


(

)

2

- (1 + Da)2

(21)

]

(22)

)

xe - ye (1 + Da) xe - xaf

xe - xaf xe - ye

)

2

- (1 + Da)

2

The coupled reactor separator system can be destabilized by delay when Da lies in the range of 0 < Da < [(xaf - 2xe + ye)/(xe - ye)]. Hence, the delay can destabilize the coupled system for a fixed vapor and liquid composition in the flash when xe > ye and xe < xaf. The parameters considered to investigate the dynamic behavior of the coupled reactor separator system are ye )

0.2, xe ) 0.4, xaf ) 0.9, and Da ) 0.5. In the absence of delay, the coupled CSTR-flash system is always stable for all values of the bifurcation parameter Da. The dependency of mole fraction of reactant A (z) on Da obtained from DDEBIFTOOL using the aforementioned parameters is shown in Figure 2. The roots of the characteristic equation at its steady state, and the real parts of the roots of the characteristic equation for the bifurcation parameter Da, is depicted in Figure 3a and Figure 3b, respectively. Here, when τ ) 1, the coupled system exhibits stable solution and is shown in Figure 3a and Figure 3b. As delay is increased from 1 to 2, the system becomes dynamically unstable and is shown in Figure 3c and Figure 3d. Here, the complex conjugate pairs of eigen values cross the imaginary axis and are depicted in Figure 3c. The increased dynamic instability in the system is solely due to the composition delay between the reactor and the separator. In this control strategy, delay has a significant effect on the dynamic behavior of the coupled system, as a result of the recycle flow rate L from the flash to the CSTR being less than the effluent flow rate F leaving from the reactor. The real parts of the roots of the characteristic equation on Da are shown in Figure 4 for a delay in the range of 3-8. The coupled system has only one region of dynamic instability for Da when τ ) 2-4 (see Figures 4a and 4b). When τ ) 5, another new region of stability branch is formed, along with the first branch, as shown in Figure 4c. When the delay lies between 5 and 7,

Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008 3927

Figure 9. Dependency of dimensionless critical delay on r (left) and frequency of oscillation versus dimensionless critical delay (right) for the variable reactor molar holdup case (ye ) 0.2, xe ) 0.4, xaf ) 0.9, Da ) 0.5, and r ) 1.5): (a and b) τ ) 2-11, (c and d) τ ) 12. Solid line represents data obtained using DDE-BIFTOOL, and symbol denote data obtained using eqs 24–27. I ) first branch and II ) second branch.

the system has two regions of dynamic instability, as shown in Figure 4c-e, and the first stability branch is fully unstable. At τ ) 8, the coupled system exhibits three regions of dynamic instability, as shown in Figure 4f. Furthermore, any increase in the delay increases the regions of dynamic instability. The dependency of critical delay and frequency of oscillation on the values of Da is shown in Figure 5. When 2 e τ e 4, the system has a single critical delay boundary and is shown in Figure 5a. Here, the coupled system becomes dynamically unstable when the chosen delay is greater than the critical delay for a fixed Da. The coupled system is stable when the chosen delay is less than the critical delay. Similarly, when 5 e τ e 7, the coupled system has two regions of critical delay boundaries, as shown in Figure 5b, in which the coupled system becomes unstable when the delay is higher than the critical delay of the first boundary. At τ ) 8, the coupled system exhibits three regions of critical delay boundaries (Figure 5c) for the values of Da. The dependency of frequency of oscillation with critical delay for the three possible cases is shown in Figures 5d-f. In this case, the frequency of oscillation is always a positive value for the parameters considered. The frequency of oscillation and critical delay obtained from the DDE-BIFTOOL package show good agreement with the same parameters obtained using eqs 21 and 22. The dynamic instability of the delay system is verified numerically by integrating the delay differential equation using

the Solver Toolbox dde23, which is available in Matlab.17 The dynamic evolution to a state of sustained oscillations of the system for Da ) 0.5, with τ ) 2.5 and τ ) 6.0, is shown in Figures 6a and 6b. In this work, it is demonstrated that the delay can destabilize the coupled system when the proposed control strategy of Luyben1,18 is implemented. 4.3. Fixed F0, F Case with Variable Molar Holdup. The dynamic instability of the system is determined for the parameters ye ) 0.2, xe ) 0.4,xaf ) 0.9, Da ) 0.5, and r ) 1.5. In the absence of delay, the coupled system is always stable for all values of r. The dependency of the mole fraction of reactant A (z) on r obtained from DDE-BIFTOOL is shown in Figure 7. The recycle flow rate (L) from the flash to the CSTR is given by z(t ′ -τ ′ ) - ye (23) L)F xe - ye

(

)

Here, the ratio of the recycle flow rate to the reactor effluent flow rate (L/F) lies in the interval of 0 to 1, when ye < z < xe. As discussed previously, the feasible operation of the coupled system is possible when L < F. The real parts of the roots of the characteristic equation versus r plots obtained from DDEBIFTOOL are shown in Figure 8. When τ ) 1, the coupled system is stable and is shown in Figure 8a. For τ ) 2-12, the coupled system exhibits delay-dependent stability and is shown in Figures 8b-d. Here, the first stability branch is fully unstable when τ > 4.

3928 Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008

of the system obtained using DDE-BIFTOOL is similar to the same values that were obtained by solving eqs 24–27. The evolution of the mole fraction of reactant A in the reactor to an oscillatory state for r ) 1.5 and τ ) 5.0 is shown in Figure 10. 5. Conclusions

Figure 10. Mole fraction of reactant A (z) versus time (t) plot, showing the evolution to sustained oscillations for the variable reactor molar holdup case (ye ) 0.2, xe ) 0.4, xaf ) 0.9, Da ) 0.5, r ) 1.5, and τ ) 5).

The dependency of the critical delay τ on r is shown in Figure 9. The system becomes unstable as the delay becomes greater than the critical delay boundary. When 2 e τ e 11, the system has only one critical delay boundary, as shown in Figure 9a. For τ ) 12, the coupled system exhibits another critical delay boundary, along with the first critical delay boundary, and is depicted in Figure 9c. As discussed previously, the coupled system exhibits dynamic instability as a result of composition delay between the reactor and the separator. Here also, delay has a significant effect on the dynamics of the coupled system, as a result of L < F. The dependency of the frequency of oscillation on the critical delay is shown in Figures 9b and 9d. Alternatively, the critical delay and the frequency of oscillation are determined by solving the steady-state equations of the mole fraction of reactant A (z) and the reactor molar holdup MR, along with the roots of the characteristic equation of the multivariable delay differential equations:

( ) ( ) ( )( ) [( ) ( )( ) ] ( ){ ( ) } [ ] ( ){ } 0)1+r

0)

xaf

M/R

+

z - ye r xe / x MR e - ye

z - xe xe - ye

1 + M/R

z - ye r + Da z / MR xe - ye

0 ) -ω2 -

r ω sin(ωτ) + / x ye MR e

(25)

1 [(xaf - xe) + M/2 R

r(xe - z)]

0)ω

(24)

1 ( r r - cos(ωτ)) + Da xe - ye M/R

cos(ωτ) (26)

1 [(xaf - xe) + M/2 R

r(xe - z)]

sin(ωτ) (27)

Equations 24 and 25 are the steady-state equations obtained from eqs 8 and 9 for the molar holdup MR of the reactor and the mole fraction of reactant A, respectively. The roots of the characteristic equation are determined by substituting a value of m ) 1 into eq 16b. The real and the imaginary parts of the roots of the characteristic equation are given in eqs 26 and 27, respectively. The critical delay and the frequency of oscillation

The effect of delay on the stability of a coupled continuously stirred tank reactor (CSTR)-flash system is examined in this work, using idealized models. The dynamic behavior of the coupled system operated under three different control strategies is governed by the delay differential equation. An elementary first-order isothermal irreversible reaction is considered to study the control strategies under which delay can destabilize the steady state of the system using the DDEBIFTOOL package. It is shown that the system exhibits delay-independent stability when the fresh feed flow rate (F0) is flow-controlled. The coupled system is destabilized by delay when the reactor effluent flow rate (F) is flow-controlled or the molar holdup MR in the reactor is allowed to vary. The critical delay boundary and the frequency of oscillation obtained from the DDE-BIFTOOL package are similar to the same values that have been obtained from analytical expressions proposed in the literature. Nomenclature A ) reactant B ) product Da ) Damko¨hler number F0 ) fresh feed flow rate entering into the reactor (mol/s) F ) effluent flow rate from the reactor (mol/s) I ) identity matrix J ) Jacobian with respect to the instantaneous state variable Jτ ) Jacobian with respect to the delayed state variable k ) kinetic constant (s-1) L ) recycle flow rate from the reactor to the separator (mol/s) MR ) reactor molar holdup; MR ) FmVR (mol) MR* ) dimensionless reactor molar holdup Mch ) characteristic molar holdup r ) bifurcation parameter; r ) F/F0 t′ ) space time (s) t ) dimensionless space time V ) vapor flow rate from the flash (mol/s) xaf ) mole fraction of reactant A in the fresh feed xe ) mole fraction of reactant A in the recycle stream ye ) mole fraction of reactant A in the distillate stream from the flash z ) mole fraction of reactant A in the reactor τ ) dimensionless critical delay τ′) critical delay (s) ω ) frequency of oscillation (rad/s) ψ ) parameter λ ) roots of the characteristic equation

Literature Cited (1) Luyben, W. L. Snowball effects in reactor separator process with recycle. Ind. Eng. Chem. Res. 1994, 33, 299. (2) Pushpavanam, S.; Kienle, A. Nonlinear behavior of an ideal reactor separator network with mass recycle. Chem. Eng. Sci. 2001, 56, 2837. (3) Kiss, A. A.; Bildea, C. S.; Dimian, A. C.; Iedema, P. D. State multiplicity in CSTR-separator-recycle polymerization systems. Chem. Eng. Sci. 2002, 57, 535. (4) Wu, K. L.; Yu, C. C. Reactor/Separator process with recycle;1. Candidate control structure for operability. Comput. Chem. Eng. 1996, 20, 1291.

Ind. Eng. Chem. Res., Vol. 47, No. 11, 2008 3929 (5) Kumar, A.; Daoutidis, P. Nonlinear dynamics and control of process systems with recycle. J. Process Control 2002, 12, 475. (6) Bildea, C. S.; Dimian, A. C. Fixing flow rates in recycle systems: Luyben’s rule revisted. Ind. Eng. Chem. Res. 2003, 42, 4578. (7) Seki, H.; Naka, Y. A. Hierarchical controller design for a reactor/ separator system with recycle. Ind. Eng. Chem. Res. 2006, 45, 6518. (8) Kiss, A. A.; Bildea, C. S.; Dimian, A. C. Design and control of recycle systems by nonlinear analysis. Comput. Chem. Eng. 2007, 31, 601. (9) Denn, M. M.; Lavie, R. Dynamics of plants with recycle. Chem. Eng. J. 1982, 24, 55. (10) 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. (11) Engelborghs, K.; Luzyanina, T.; Roose, D. Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL. ACM, Trans. Math. Soft. 2002, 28, 1. (12) Murray, J. D. Mathematical Biology; Springer Verlag: Berlin,1993. (13) Epstein, I. R.; Pojman, J. A. An Introduction to Non-linear Chemical Dynamics, Oscillations, WaVes, Patterns and Chaos. Oxford University Press: New York, 1998.

(14) Engelborghs, K.; Luzyanina, T.; Roose, D. Numerical bifurcation analysis of delay differential equations. J. Comput. Appl. Math. 2000, 125, 265. (15) Engelborghs, K. DDE-BIFTOOL: A Matlab package for bifurcation analysis of delay differential equations. Technical Report TW-305 Department of Computer Science, K.U. Leuven, Leuven, Belgium, 2000. (16) Engelborghs, K.; Luzyanina, T.; Samaey, G. DDE- BIFTOOL v.2.00 User Manual: A Matlab Package for Bifurcation Analysis of Delay Differential Equations. Technical Report TW-330 Department of Computer Science, K.U. Leuven, Leuven, Belgium, 2001. (17) Shampine, L. F.; Thompson, S. SolVing DDE with Matlab; Mathematics Department, Southern Methodist University: Dallas, TX, 75275. (18) Luyben, W. L.; Luyben, M. L. Essentials of Plant Control; McGraw-Hill, International Edition: Singapore, 1997.

ReceiVed for reView July 25, 2007 ReVised manuscript receiVed January 8, 2008 Accepted February 29, 2008 IE071010C