Bifurcation and Stability Analysis of an Ethylene Oxide Reactor System

The stable control regions of the ethylene oxide reactor system were developed as a function of operating ... Control of Rotating Wave Trains in a Loo...
0 downloads 0 Views 192KB Size
Download PDF
Ind. Eng. Chem. Res. 2003, 42, 3285-3293

3285

Bifurcation and Stability Analysis of an Ethylene Oxide Reactor System Kishor G. Gudekar and James B. Riggs* Department of Chemical Engineering, Texas Tech University, Lubbock, Texas 79409-3121

In this work, an open-loop and closed-loop nonlinear stability analysis of an industrial ethylene oxide reactor is performed. The simulator, which consists of a gas-gas heat exchanger, a multitubular fixed-bed reactor, a steam generator, and a separation system, was benchmarked with data from an industrial ethylene oxide reactor system. A two-dimensional heterogeneous model was developed for the multitubular fixed-bed reactor. The steady-state operability problem (open-loop bifurcation) was addressed by using nonlinear bifurcation techniques. The stable control regions of the ethylene oxide reactor system were developed as a function of operating temperature, catalyst activity, controller tuning, and disturbance direction and magnitude. A reactor runaway event was simulated, from which it is shown that eliminating the oxygen in the feed to the reactor can be used to prevent reactor runaway. Introduction It is well-known that the input/output multiplicities of chemical reactors can have an important effect on the operational difficulties of such processes.1 Bifurcation theory has been recognized as a very useful tool for addressing the nonlinear behavior of processing systems subject to the variation of a parameter.2 Using such bifurcation methods, it becomes feasible to numerically detect input multiplicities, output multiplicities, isolas, limit cycles, etc. Jensen et al.3 applied static and Hopf bifurcation theory to partial differential equations (PDEs) for the special case of a first-order, irreversible reaction in a tubular reactor with axial dispersion and classified and summarized the bifurcation behavior in parameter space plots. Chang4 presented an analysis of the various types of bifurcation that are caused by a conventional, singleinput-single-output PID controller on an isothermal bioreactor. The steady-state multiplicity behavior of a nonisothermal heterogeneous reactor model with axial dispersion was examined by Juncu et al.5 This model took into account both external and internal heat and mass transfer with composition and temperature variations in the axial and radial directions for the fluid and solid phases. An irreversible first-order reaction in the adiabatic, nonisothermal, and nonadiabatic regimes was considered. Bifurcation theory and numerical continuation techniques were used by Wagialla et al.6 to investigate the complex static and dynamic characteristics of fixed-bed reactors modeled by a heterogeneous cell model. A single exothermic first-order reaction was considered. Varma et al.7 generated a reactor operation diagram for ethylene oxide process consisting of the runaway region along with an optimal curve that gives the optimal reactor outlet yield of ethylene oxide. Baratti et al.8 theoretically studied ethylene epoxidation over silver catalyst in a heterogeneous, nonadiabatic, nonisothermal plug-flow reactor and showed that reactor runaway can be avoided using Dirac catalysts as compared to a uniform catalyst distribution. * To whom correspondence should be addressed. Tel.: 806742-1765. Fax: 806-742-1765. E-mail: [email protected].

Quina et al.9 presented a steady-state analysis of the region of parametric sensitivity and the range of operating conditions leading to temperature runaway for a fixed-bed reactor (the selective oxidation of methanol to formaldehyde) where the catalytic bed is partially diluted with inert packing. Garcia et al.10 studied the steady-state nonlinear bifurcation behavior of a continuous stirred-tank reactor that produced high-impact commercial polystyrene. A complete bifurcation analysis of a general steady-state two-dimensional catalytic monolith reactor model that accounted for temperature and concentration gradients in both the axial and radial directions was performed by Balakotaiah et al.11 Thus, most of the studies in this area are centered on the steady-state bifurcation analysis of reactors. Such steadystate bifurcation analyses aid in understanding the input/output multiplicity in the reactor and the stable/ unstable operating regions with respect to certain physicochemical parameters. For a fixed-bed reactor to operate in a reliable and safe manner, not only is a steady-state nonlinear bifurcation analysis important, but a closed-loop stability analysis is also important. The aim of this work is to provide a first look into the operability problems encountered with ethylene oxide reactors and to perform an open-loop and closedloop stability analysis using a model benchmarked against an industrial ethylene oxide reactor system. Ethylene oxide is an important petrochemical intermediate and is used for the production of glycol and poly(ethylene glycol). The operation of an ethylene oxide reactor is critical because the reactor can generate 11 times as much heat in a runaway condition as under normal operating conditions. Therefore, safety issues for ethylene oxide reactor systems are critical as industry also tries to operate reactors in an economically advantageous manner. That is, a higher reactor operating temperature will generally result in a greater profit for the process, but it will also increase the tendency of a reactor to exhibit thermal runaway. Therefore, the stable temperature control of an ethylene oxide reactor is economically important. In this work, an analysis of the stable control region of the system is developed as a function of the operating temperature, catalyst activ-

10.1021/ie020874p CCC: $25.00 © 2003 American Chemical Society Published on Web 06/04/2003

3286

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003

Figure 1. Schematic of the ethylene oxide process.

ity, controller tuning, and disturbance direction and magnitude. Elementary catastrophe theory can be used to detect analytical conditions under which input/output multiplicities could emerge for a given reactor system. However, catastrophe theory requires that the entire mathematical model be collapsed into a single algebraic equation. Therefore, such an approach is not practical for large-scale models. Because of the complexity of the ethylene oxide process model equations, a purely numerical procedure is used to characterize the multiplicity behavior of the ethylene oxide process. A continuation algorithm is used to generate a branch of steadystate solutions with respect to a continuation parameter to obtain bifurcation diagrams. In principle, this technique enables us to start from a known solution and continuously compute solutions along a chosen branch.12

The following equations were obtained for the production rates of ethylene oxide and carbon dioxide

Process Description Figure 1 shows a schematic of the ethylene oxide process considered here. The apparatus studied consists of a gas-gas heat exchanger, a multitubular fixed-bed reactor, a steam generator, and a separation system. The exothermic heat of reaction from the reactor is removed by passing coolant on the shell side of the reactor tubes. A portion of the heated coolant is passed through a steam generator to produce steam, and the total coolant stream is recycled back to the shell side of the reactor. A single-loop PID control system uses the flow rate of the coolant that is passed through the steam generator to maintain the inlet temperature of the coolant to the reactor. The interesting features of this system are (1) a nonlinear PDE description, (2) a runaway reaction that produces carbon dioxide, and (3) a tradeoff between selectivity and conversion.

Mathematical Modeling (i) Fixed-Bed Reactor Modeling. A two-dimensional heterogeneous dynamic model was developed for the ethylene oxide reactor that is based on the following assumptions: (1) Axial dispersion is negligible. (2) Pressure drop is negligible. (3) No concentration gradient exists between the catalyst particle and the gas. (4) The velocity profile is flat. (5) No intraparticle resistance is present (i.e., it is assumed that the effectiveness factor for each of the reactions is unity). Each of the above assumptions was verified, and the details of the verification process can be found in Gudekar.14 The following dimensionless variables were used to obtain the dimensionless form of the reactor equations

Kinetics Borman and Westerterp13 studied the kinetics of the selective oxidation of ethylene over a silver R-alumina catalyst. The relevant reactions are the partial oxidation of ethylene to give ethylene oxide and the complete oxidation of ethylene to produce carbon dioxide and water

1 C2H4 + O2 f C2H4O 2 C2H4 + 3O2 f 2CO2 + 2H2O

Rj ) kjrPC2H4PO2nj j j 1 + KCj 2H4PC2H4 + KCO P + KH P + KjEOPEO 2 CO2 2O H2O

j ) 1, 2 The rate constants, kjr, are of Arrenhieus form; the j , and KjEO, are constants, adsorption terms, KCj 2H4, KCO 2 j but KH2O is modeled as an exponential function of temperature. For numerical values of the parameters used in this study, the reader is referred to Borman and Westerterp.13

C2H4* ) CO2* )

C2 H4 0

C2H4

CO2 0

C2H4

T* )

, O2 * )

, H2O* )

O2 0

C2H4

H2O C2H4

0

, EO* )

, CH4* )

(T - TR) (Ts - TR) , T/s ) 10TR 10TR

(Tc - TR) (Tm - TR) , T/m ) 10TR 10TR z u r r* ) , z* ) , τ ) t R L L

T/c )

EO C2H40 CH4 C2H40

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3287

Using the above dimensionless variables, the following system equations can be obtained

Material Balances

(

)

∂C2H4* ∂C2H4* DrL ∂ ∂C2H4* )+ 2 r* ∂τ ∂z ∂r* R ur* ∂r* FbL (r1 + r2) (1) uC2H40

(

)

∂O2* ∂O2* DrL ∂ ∂O2 )+ 2 r* ∂τ ∂z* ∂r* ∂r* R ur* FbL (0.5r1 + 1.5r2) (2) uC2H40 DrL ∂ FbL ∂EO* ∂EO* ∂EO* )+ 2 r* + r1 ∂τ ∂z* ∂r* ∂r* R ur* uC2H40 (3)

(

∂CO2* ∂τ

)-

)

(

)

DrL ∂ ∂CO2 FbL ∂CO2* + 2 r* + (2r2) ∂z* ∂r* R ur* ∂r* uC2H40 (4)

(

)

∂H2O ∂H2O* ∂H2O* DrL ∂ r* )+ 2 * + ∂τ ∂z* ∂r* R ur ∂r* FbL

(

(2r2) (5) uC2H40

)

∂CH4* ∂CH4* ∂CH4* DrL ∂ r* )+ 2 ∂τ ∂z* ∂r* R ur* ∂r*

(6)

)

4htL ∂T* ∂T* )+ (T/ - T*) + ∂τ ∂z* FgCpgdiu c hgsaL λrL ∂ ∂T* (T* - T/s ) (7) r* 2 ∂r* FgCpgu R F C ur* ∂r* g

pg

( )

Energy Balance for a Catalyst Particle ∆H1L ∂T/s hgsaL (T* - T/s ) + r + ) ∂τ FbCpsu(1 - ) 10CpsuTR(1 - ) 1 ∆H2L r (8) 10CpsuTR(1 - ) 2 Energy Balance for the Metal Tubes dT/m 4dohsL (T/c - T/m) ) 2 2 dτ CpmFpmu(do - di ) 4dihiL 2

at r* ) 0

dC/j dT* ) 0, )0 dr* dr*

at r* ) 1

dC/j dT* ) 0, -λw ) hwR(T/m - T/c ) dr* dr*

at t ) 0

z* ) 0, T* ) T/0, T/s ) T/s0, / T/m ) T/m0, T/c ) T/c0, C/j ) Cj0

The partial differential equations describing the fixedbed reactor were converted into ordinary differential equations by orthogonal collocation.15 The radial temperature derivatives were approximated by a quadratic polynomial (i.e., one internal collocation point), and the radial concentration derivatives were approximated by a quartic polynomial (i.e., two internal collocation points). A cubic polynomial (i.e., three internal collocation points) was used for the approximation of the axial derivatives. Six subintervals in the axial direction were used to model the entire reactor. The reduced dynamic ordinary differential equations were stiff, and therefore, they were solved using a stiff LSODE integrator.16 The developed model predicts composition and temperature profiles in the axial and radial directions of the reactor as a function of time. (ii) Model Equations for the Steam Generator. The steam generator, which is a shell-and-tube heat exchanger, is modeled as a distributed parameter system. Following are the two energy balance equations for the steam generator

Dimensionless Energy Balance for the Metal Tubes

(Note that CH4 is used as an inert component in this system, and its concentration in the reactor does affect the partial pressures of the reactants.)

(

the following initial and boundary conditions

dT ˆ /m 4dˆ ohˆ sL ˆ /m) (T ˆ /stm - T ) dτ C ˆ pmFˆ pmu(dˆ o2 - dˆ i2) 4dˆ ihˆ iL ˆ /c ) (11) (T ˆ /m - T C ˆ pmFˆ pmu(dˆ o2 - dˆ i2) Dimensionless Energy Balance for the Shell-Side Coolant ∂T ˆ /c uˆ cL ∂T ˆ /c 4hˆ snˆ dˆ oL ˆ /c ) (12) (T ˆ /m - T )+ ∂τ uL ˆ ∂z* C ˆ pcFˆ cu(dˆ r - nˆ dˆ o2) The above dimensionless equations were obtained using the following dimensionless variables

T ˆ /stm ) T ˆ /m )

2

(T/m - T*) (9)

CpmFpmu(do - di ) Shell-Side Energy Balance

uc ∂T/c 4hsndoL ∂T/c (T/m - T/c ) )+ ∂τ u ∂z CpcFcu(dr2 - ndo2)

(10)

The above partial differential equations are subject to

(T ˆ stm - TR) (T ˆ c - TR) , T ˆ /c ) 10TR 10TR

(T ˆ m - TR) z u , z* ) , τ ) t 10TR L L ˆ

The partial differential equations describing the steam generator were converted into ordinary differential equations by the application of the method of lines.17 The axial temperature derivatives were approximated by finite differences. The set of dynamic ordinary differential equations was solved using an LSODE integrator.16 (iii) Model Equations for the Gas-Gas Heat Exchanger. For the modeling of the gas-gas heat exchanger, the effectiveness method18 was used to

3288

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003

calculate the steady-state exchanger exit temperatures, and then these temperatures were lagged by a firstorder time constant (6 s) to model the process dynamics. (iv) Model Equations for the Separation System. In the separation system, ethylene oxide and carbon dioxide are removed, and the unreacted ethylene and oxygen are recycled back to the fixed-bed reactor to be mixed with fresh ethylene, oxygen and inert methane. The resulting mixture is then passed through the gasgas heat exchanger, where it is heated to the reactor temperature. The separation system was modeled as a first-order process as shown below. The value of the effective time constant,τsp (3 min), was selected on the basis of industrial experience.

dC2H4 (C2H4ss - C2H4) ) dt τsp dO2 (O2ss - O2) ) dt τsp ss dEO (EO - EO) ) dt τsp

dCO2 (CO2ss - CO2) ) dt τsp dT3 (Tss 3 - T3) ) dt τsp (v) Steady-State Model Benchmarking. The steadystate version of each of the above models was benchmarked against industrial data by minimizing the weighted error between the model-predicted values (i.e., compositions and temperatures) and the industrial data. The industrial data came from a period during which fresh catalyst was used; therefore, catalyst deactivation did not require consideration. The model predicted the industrial data quite well, with an average relative error of 0.87% for composition predictions and an average absolute error of 0.3 °C for temperature predictions. (vi) Catalyst Deactivation Model. For ethylene oxide reactors, it is believed that the deactivation occurs as a result of sintering of the catalyst.19 Therefore, the catalyst activity can be expressed as a nonlinear function of the catalyst age, tyr, and the inlet shell-side temperature Tin c

a ) exp(-Rtyr)

for tyr e tcon

a ) exp(-Rtyr) exp(-βTin c )

for tyr > tcon, Tin c < Tcon

a ) exp(-Rtyr) exp(-γTin c )

for tyr > tcon, Tin c g Tcon

(vii) Catalyst Deactivation Model Benchmarking. The catalyst deactivation model parameters (R, β, γ) were identified by minimizing the weighted error between the model-predicted values (i.e., composition and temperatures) and the industrial data over the life of a batch of catalyst. The catalyst deactivation model predicted the industrial data quite well over the life of the catalyst, with an average relative error of 2.0% for composition predictions and an average absolute error of 0.65 °C for temperature predictions.

Figure 2. Bifurcation diagram using the flow through the steam generator as the continuation parameter.

Open-Loop Nonlinear Bifurcation Analysis (i) Effect of Manipulated Variable. Figure 2 shows the bifurcation diagram for the reactor coolant inlet temperature, maximum reactor temperature, and ethylene conversion (%) using the coolant flow through the steam generator as the continuation parameter. A continuation algorithm was used to generate the continuation (bifurcation) diagram (i.e., for each value of the coolant flow through the steam generator, a steadystate solution was found). Initially, the bifurcation diagram was obtained by solving the nonlinear algebraic steady-state model equations for a specified coolant flow through the steam generator. The steady-state solution thus found was used as an initial guess for the next computation. For this system of equations to converge, a number of initial guesses are required (more than 5). Therefore, to reduce the computational time, the problem was reformulated to obtain the bifurcation diagram. An additional equation for the reactor inlet coolant temperature was added to the existing nonlinear equations. This equation describes the difference between the reactor inlet coolant temperature and the specified reactor inlet coolant temperature as shown below in in f(Tin c ) ) Tc - (Tc )sp

Rather than specifying the coolant flow through the steam generator, the reactor coolant inlet temperature (which is a set point for the reactor inlet coolant temperature control system) was specified, and the resulting nonlinear algebraic equations were solved using MINPACK.20 This formulation of the problem was found to be more efficient than the previous one. Using the previous method, only one branch of solutions (stable or unstable) was obtained. To find the other branch of solutions, a different initial guess was required. In contrast, using the latter formulation method starting from one initial guess, both the branches of solutions (stable and unstable) were always obtained. Under nominal operating conditions the reactor displays output multiplicities, i.e., for a given coolant flow through the steam generator, there are two different reactor operating conditions. The nominal upper steadystate temperature is unstable, while the lower steadystate temperature is stable, indicating saddle node bifurcation. In saddle node bifurcation, there is only a

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3289

Figure 4. Temperature of the catalyst particle and its normalized sensitivity with respect to the flow through the steam generator. Figure 3. Bifurcation diagram using the mole fraction of CO2 in the reactor feed as the continuation parameter.

single turning point, and this turning point marks the boundary between stable and unstable behavior. One would expect an odd number of steady-state solutions to exist for this reactor system. Because two steady-state solutions were found, a third should exist. In this case, there were only two were found, and the third solution is expected to be outside the industrially relevant operating region because a continuation variable (e.g., the fractional flow through the steam generator) was varied over its full range to locate the steadystate operating points. (ii) Effect of Disturbance. The open-loop bifurcation behavior of the ethylene oxide reactor process was analyzed with respect to process disturbances. A change in the carbon dioxide concentration in the reactor feed was found to be the most severe disturbance to the ethylene oxide reactor inlet coolant temperature control loop. The same procedure that was described previously was used to obtain the bifurcation diagram based on the inlet carbon dioxide composition as the continuation parameter (Figure 3). As in the effect of the manipulated variable, the reactor displays output multiplicities (saddle node bifurcation) with respect to the mole fraction of CO2 in the reactor feed. Thus, the bifurcation diagrams (Figures 2 and 3) can be used to understand changes in stability on a given branch of solutions. The upper branch, which is unstable, can be stabilized by feedback control. All of the previous bifurcation diagrams were developed for fresh catalyst (catalyst activity is unity). As the catalyst activity decreases with time, the operating temperature is increased to compensate for the loss of activity. In the next section, an open-loop runaway boundary is developed with respect to catalyst activity. To protect the confidentiality of the industrial EO operator, we are not able to provide actual values of the operating conditions (e.g., temperature and composition) in the figures and tables. (iii) Runaway Boundary. The magnitude of the hotspot temperature in the reactor must be bounded within specific limits because it might seriously affect reactor safety and performance. The hot-spot temperature depends on the process operating conditions (e.g., feed flow, temperature, and composition) and process parameters (e.g., reaction kinetic parameters, viscosity of the gas, and thermal conductivity of the catalyst). For specific values of the system parameters, the hot spot might experience high sensitivity to one or more of the operating conditions or system parameters. In this case, the reactor is said to operate in the runaway or

parametrically sensitive region. In practical applications, it is desirable to avoid this operating region for safety of the process. This provides the motivation to develop a runaway boundary for the ethylene oxide reactor. In particular, we will identify the runaway region for this reactor by applying the generalized runaway criterion21 using the catalyst temperature at the exit of the reactor, θ/p, as the objective. For this, we need to define the objective sensitivity, s(θ/p;φ)

s(θ/p;φ)

dθ/p ) dφ

where φ represents the model input parameter or operating condition such as coolant flow through steam generator. However, the more appropriate quantity in sensitivity analysis is the normalized objective sensitivity, S(θ/p;φ). The normalized objective sensitivity of the catalyst temperature maximum, θ/p, along the reactor length is given by /

S(θ/p;φ) )

φ φ dθp s(θ/p;φ) ) / / θp θp dφ

In the present ethylene oxide process, the axial temperature profile is monotonically increasing, and so, the maximum catalyst temperature value considered in the sensitivity analysis is that at the reactor outlet. The critical conditions for reactor runaway, according to the generalized sensitivity criterion, are then identified as the conditions for which the magnitude of the normalized objective sensitivity, S(θ/p;φ), is maximized.21 As shown in Figure 4, the magnitude of the normalized sensitivity S(θ/p;φ) for the maximum reactor temperature increases with decreasing flow through the steam generator, and the magnitude of the normalized sensitivity is largest at the bifurcation point. Thus, at the bifurcation point, the catalyst temperature is very sensitive to the change in the coolant flow through steam generator. The runaway boundary is coincident with the bifurcation point of the multiple-steady-state regions. Note that the sign of the normalized sensitivity value has a particular meaning. A positive (negative) value of the normalized sensitivity of the maximum temperature with respect to the coolant flow through the steam generator indicates that the maximum temperature increases (decreases) as the magnitude of coolant flow increases. Thus, in this case, the normalized sensitivity is negative, which indicates that the transition from

3290

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003

Figure 7. Closed-loop stability region for positive and negative changes in the carbon dioxide composition (mol %) in the reactor feed. Figure 5. Locus of bifurcation points for different catalyst activities.

Figure 8. Effect of the detuning factor on the closed-loop stability boundary. Figure 6. Boundary of the open-loop runaway region.

nonrunaway to runaway behavior occurs as the coolant flow through the steam generator is decreased. Figure 5 shows the effect of catalyst activity on the locus of bifurcation points. Note that the bifurcation point marks the boundary between stable and runaway conditions. Figure 6 shows the open-loop runaway region as a function of catalyst activity. From Figure 6, one can see that, as the catalyst activity decreases, the temperature at which bifurcation occurs increases. Any arbitrary size disturbance while the process is operating at the runaway boundary will cause the reactor to runaway under open-loop conditions. Closed-Loop Nonlinear Stability Analysis The ethylene oxide reactor is normally operated at the higher coolant inlet temperatures (open-loop unstable operating point) because it results in a higher ethylene conversion than would be obtained at a lower coolant inlet temperature (open-loop stable operating point). Although the operating point is open-loop unstable in implementation, a feedback controller stabilizes it. Following is an analysis of the effects of various factors on closed-loop stability. (i) Effect of Disturbances. A change in the inlet carbon dioxide composition for the reactor feed is the most severe disturbance for the reactor inlet coolant temperature control loop because carbon dioxide affects the split between partial and complete combustion of ethylene. The magnitude of the inlet carbon dioxide composition disturbance was changed in the positive and negative directions for different values of the control-loop dead time to determine the maximum disturbance that the PI controller was able to handle. The control-loop dead time (θ1 + θ2), as shown in Figure 1, is the transportation delay between the control valve

on the coolant flow and the location where the reactor inlet coolant temperature is measured. Figure 7 shows the runaway region for positive and negative changes in the carbon dioxide reactor inlet composition. For a positive change in the carbon dioxide reactor inlet composition with a dead time of less than 20 s, the system becomes unstable for almost a +10 mol % absolute change in the carbon dioxide inlet composition. This occurs because the increased carbon dioxide partial pressure (which is due to increased inlet carbon dioxide mole percentage) favors the partial oxidation of ethylene and does not favor the complete oxidation of ethylene, which has a large heat of reaction. For a negative change in the carbon dioxide reactor inlet composition, the system is more prone to instability because the decrease in the partial pressure of the carbon dioxide in the reactor favors the complete oxidation reaction, which has a much larger heat of reaction. From Figure 7, we conclude that a negative change in the carbon dioxide reactor inlet composition is a much more severe disturbance than a positive change. Therefore, for the remaining studies, only the negative change in the carbon dioxide reactor inlet composition is considered as a disturbance. (ii) Controller Tuning. Here, we study how the aggressiveness of the controller affects the runaway region. First, the controller is tuned for set-point changes using a 1/6 decay ratio as the tuning criterion by adjusting the detuning factor.22 The resulting detuning factor was typically in the range of 1.0-1.3. Then, the detuning factor was varied from 0.5 to 5.0, and each time, the magnitude of the disturbance was varied until controller instability was attained. Note that, as the detuning factor is decreased, the controller aggressiveness is increased. Figure 8 shows the effect of the detuning factor on the stability region for a dead time of 8 s in the control loop. The system is shown to

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3291

Figure 11. Comparison between nominal industrial operation and the optimal reactor inlet coolant temperature profile. Figure 9. Effect of catalyst activity on the stability boundary.

Figure 12. Response of the outlet gas temperature when runaway is observed.

Figure 10. Effect of the reactor inlet coolant temperature on the closed-loop stability region.

be more sensitive to reactor runaway as the coolant temperature controller is detuned. (iii) Effect of Catalyst Activity. Figure 9 shows the effect of the reactor inlet coolant temperature on the runaway boundary for different catalyst activities. These results are based on maintaining a fixed ethylene oxide production rate; therefore, the operating temperature increases as the catalyst activity decreases. The previously described procedure was applied to obtain the runaway boundaries for a range of catalyst activities. (iv) Effect of Reactor Inlet Coolant Temperature. The following procedure was used to identify the runaway region as a function of reactor inlet coolant temperature. First, for a specified reactor inlet coolant temperature, the PI controller was tuned for an inlet coolant temperature set point using a 1/6 decay ratio as a tuning criterion. Then, the disturbance magnitude was varied until the controller became unstable. From Figure 10, it can be seen that, at higher reactor inlet coolant temperatures, the control system becomes unstable for relatively small decreases in the CO2 feed composition. Modeling of Reactor Runaway The operating temperature of an ethylene oxide reactor is important because, at higher operating temperatures, the selectivity for ethylene oxide decreases, and at lower temperatures, the conversion of ethylene to ethylene oxide decreases. Therefore, there is an optimal operating temperature that provides the optimal economic tradeoff between conversion and selectivity. The optimal temperature profile shown in Figure 11 was obtained by maximizing the process profit

function over the entire catalyst life for a fixed operating period and fixed reactor inlet conditions. The optimization study was based on a fixed feed composition and feed flow rate with the economic objective function defined as the gross profit based on the value of the products, the cost of the feed material, the cost of utilities, the cost of catalyst, and outages required to replace spent catalyst. Figure 11 also shows the normal operating temperature for the life of the catalyst. Details of the determination of the optimal temperature profile can be found in Gudekar.14 According to Figures 9 and 10, changing from the nominal operating temperature to the optimal operating temperature makes the ethylene oxide reactor more susceptible to runaway. The reactor undergoes runaway when the rate of change of heat removal with temperature is less than the rate of change of heat generation with temperature. One such runaway reactor is simulated in next section. (i) Simulated Runaway. The dynamic model that was developed for the control studies was extended to model the runaway process for an ethylene oxide process. Because the control studies considered relatively small changes in the reactor temperature, constant physical properties (e.g., gas-phase viscosity, thermal conductivities, and heat capacities) were assumed, although Arrenheius rate expressions were used. Because the reactor temperature changed by over 30 K during the modeling of a reactor runaway event, the assumption of constant physical properties was invalid. In fact, by adding the temperature dependence of the physical properties, the time to runaway was reduced by over an order of magnitude. Figure 12 shows a simulated runaway case. A disturbance (a decrease in the carbon dioxide reactor inlet concentration) was introduced after 5 s, which caused the reactor to runaway (i.e., the reactor temperature increased exponentially after a short time). The run-

3292

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003

Figure 13. Response of the manipulated variable when runaway is observed.

Figure 14. Response of the outlet gas temperature when runaway is prevented.

Figure 15. Response of the manipulated variable (flow through the steam generator) when runaway is prevented.

away was caused by the rapid increase in the complete oxidation of ethylene. Figure 13 shows the fraction of coolant flow required by the inlet coolant temperature controller for the runaway event. Note that, after 1.2 min, a constraint on the maximum flow through the steam generator was encountered. (ii) Methods of Preventing Runaway. A reactor runaway can be prevented by eliminating the make-up flow of oxygen to the reactor feed (Figure 14). The makeup oxygen is cut to zero after 1.2 min once the manipulated variable becomes saturated (Figure 15) and the reactor outlet temperature is still increasing (Figure 14). After the make-up oxygen is eliminated, the oxygen partial pressure becomes low, the reaction ceases, and the temperature drops sharply. After 1.2 min, no oxygen is present in the reactor, and the heat generated by the reaction is extinguished. Thus, elimination of the oxygen can always prevent a runaway reactor if it is applied in time.

Conclusions A nonlinear open-loop bifurcation analysis has been applied to an industrially benchmarked dynamic simulation of an ethylene oxide process indicating saddle point bifurcations. A nonlinear closed-loop stability analysis showed that the closed-loop system is most sensitive to decreases in the CO2 concentration in the feed to the reactor because such disturbances favor the complete oxidation of ethylene. Further, the nonlinear stability analysis indicates that an ethylene oxide process with a detuned temperature controller or a process operating at elevated reactor temperatures is most prone to reactor runaway. The accurate modeling of a runaway ethylene oxide reactor requires inclusion of the temperature dependences of all of the physical properties of the system. It was demonstrated that a runaway reactor can be identified when the maximum heat removal is attained and the reactor temperature continues to increase. Using simulations of a runaway reactor event, it was shown that a runaway reactor can be returned to safe and stable operation by eliminating the oxygen feed to the reactor once the reactor is identified as being in a runaway condition. Acknowledgment Support for this work from the member companies of the Texas Tech Process Control and Optimization Consortium is acknowledged. Nomenclature a ) activity of the catalyst (0 < a e 1) CH4* ) dimensionless methane concentration Cj/ 0 ) dimensionless concentration of component j at the reactor inlet CO2* ) dimensionless carbon dioxide concentration CO2ss ) amount of carbon dioxide exiting the separation system at steady state Cpc ) heat capacity of the coolant for the reactor (J/kg‚K) C ˆ pc ) heat capacity of the coolant for the steam generator (J/kg‚K) Cpg ) gas heat capacity (J/kg‚K) Cpm ) heat capacity of the metal for the reactor (J/kg‚K) C ˆ pm ) heat capacity of the metal for the steam generator (J/kg‚K) Cps ) heat capacity of the catalyst (J/kg‚K) C2H4* ) dimensionless ethylene concentration C2H40 ) ethylene concentration at the reactor inlet (mol/m3) C2H4ss ) amount of ethylene exiting the separation system at steady state di ) inside diameter of a tube in the reactor (m) dˆ i ) inside diameter of a tube in the steam generator (m) do ) outside diameter of a tube in the reactor (m) dˆ o ) outside diameter of a tube in the steam generator (m) dr ) diameter of the reactor (m) dˆ r ) diameter of the shell for the steam generator (m) Dr ) radial diffusivity (m2/s) EO* ) dimensionless ethylene oxide concentration EOss ) amount of ethylene oxide exiting the separation system at steady state hs ) shell-side heat-transfer coefficient for the reactor (W/m2‚K) hˆ s ) shell-side heat-transfer coefficient for the steam generator (W/m2‚K) ht ) tube-side heat-transfer coefficient for the reactor (W/m2‚K)

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3293 hˆ t ) tube-side heat-transfer coefficient for the steam generator (W/m2‚K) H2O* ) dimensionless water concentration Kij ) absorption rate constant for component j in reaction i (Pa-1) kir ) reaction rate constant for reaction i (mol/kg‚s‚bar) L ) length of a reactor tube (m) L ˆ ) length of a tube in the steam generator (m) n ) number of tubes nˆ ) number of tubes for the steam generator O2* ) dimensionless oxygen concentration O2ss ) amount of oxygen exiting the separation system at steady state Pj ) partial pressure of component j (Pa) s ) objective sensitivity S ) normalized sensitivity T ) bulk gas temperature for the reactor (K) T* ) dimensionless bulk gas temperature for the reactor T/0 ) dimensionless bulk gas temperature at the reactor inlet Tc ) coolant temperature for the reactor (K) T/c ) dimensionless coolant temperature for the reactor T/c0 ) dimensionless coolant temperature at the reactor inlet T ˆ /c ) dimensionless coolant temperature for the steam generator Tin c ) reactor inlet coolant temperature (K) Tcon ) constant used in the deactivation model (Tin c )sp ) reactor inlet coolant temperature set point (K) Tm ) metal temperature for a tube in the reactor (K) T/m ) dimensionless metal temperature for a tube in the reactor T/m0 ) dimensionless metal temperature at the reactor inlet T ˆ /m ) dimensionless metal temperature for the steam generator TR ) room temperature (303 K) Ts ) catalyst temperature for the reactor (K) T/s ) dimensionless catalyst temperature for the reactor T/s0 ) dimensionless catalyst temperature at the reactor inlet / T ˆ stm ) dimensionless steam temperature for the steam generator Tss 3 ) temperature of the stream exiting the separation system at steady state tyr ) catalyst life in years tcon ) constant used in the deactivation model u ) superficial tube velocity in a reactor tube (m/s) uc ) velocity of the coolant (m/s) uˆ c ) velocity of the coolant for the steam generator (m/s) z* ) dimensionless reactor length R, β, γ ) empirical constants for the catalyst deactivation model ∆Hn ) heat of reaction of the nth reaction (J/mol‚K)  ) void fraction θ/p ) maximum catalyst temperature (K) θ1, θ2 ) transport delay indicated in Figure 1 λr ) radial thermal conductivity (W/m‚K) Fb ) bulk density (kg/m3) Fc ) density of the coolant (kg/m3) Fˆ c ) density of the coolant for the steam generator (kg/m3) Fpm ) density of the metal (kg/m3) Fˆ pm ) density of the metal for the steam generator (kg/m3)

φ ) model input parameter τsp ) time constant for the separation system

Literature Cited (1) Seider, W. D.; Brengel, D. D.; Provost, A. M.; Widagdo, S. Nonlinear Analysis in Process Design. Why Overdesign to Avoid Complex Nonlinearities?. Ind. Eng. Chem. Res. 1990, 29, 805. (2) Kuznestov, Y. A. Elements of Applied Bifurcation Theory, 2nd ed.; Springer: New York, 1998. (3) Jensen, K. F.; Ray, W. H. The bifurcation behavior of tubular reactors. Chem. Eng. Sci. 1982, 37, 199-222. (4) Chang, H. C.; Chen, L. H. Bifurcation characteristics of nonlinear systems under conventional PID control. Chem. Eng. Sci. 1984, 39 (7/8), 1127-1142. (5) Juncu, G. H.; Bildea, S.; Floarea, O. Steady-state multiplicity analysis of the heterogeneous axial dispersion fixed-bed reactor. Chem. Eng. Sci. 1994, 49 (1), 123-130. (6) Wagialla, K. M.; Elnashaie, S. S. E. H. Bifurcation and complex dynamics in fixed-bed catalytic reactors. Chem. Eng. Sci. 1995, 50 (17), 2813-2832. (7) Varma, A.; Morbidelli, M.; Wu, H. Parametric Sensitivity in Chemical Systems; Cambridge University Press: Cambridge, U.K., 1999. (8) Baratti, R.; Gavriilidis, A.; Morbidelli, M., Varma, A. Optimization of a nonisothermal nonadiabatic fixed-bed reactor Using dirac type catalysts for ethylene epoxidation. Chem. Eng. Sci. 1994, 49 (12), 1925-1936. (9) Quina, M. J.; Ferreira, R. M. Q. Thermal runaway conditions of a partially diluted catalytic reactor. Ind. Eng. Chem. Res. 1999, 38, 4615-4623. (10) Verazaluce-Garcia, J. C.; Flores-Tlacuahuac, A. F.; Saldivar-Guerra, E. Steady-state nonlinear bifurcation analysis of a high-impact polystyrene continuous stirred tank reactor. Ind. Eng. Chem. Res. 2000, 39, 1972-1979. (11) Balakotaiah, V.; Gupta, N.; West, D. H. Bifurcation analysis of a two-dimensional catalytic monolith reactor model. Chem. Eng. Sci. 2001, 56, 1435-1442. (12) Kubicek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer-Verlag: New York, 1983. (13) Borman, P. C.; Westerterp, K. R. An experimental study of the selective oxidation of ethene in a wall cooled tubular packed bed reactor. Chem. Eng. Sci. 1992, 47 (9-11), 2541-2546. (14) Gudekar, K. G. Ph.D. Thesis, Texas Tech University, Lubbock, TX, 2002. (15) Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. (16) Hindmarsh, A. C. Computing and Mathematics Research Division, Lawrence Livermore National Laboratory, Livermore, CA, 1986. (17) Riggs, J. B. An Introduction to Numerical Methods for Chemical Engineers, 2nd ed.; Texas Tech University Press: Lubbock, TX, 1994. (18) Luyben, W. L.; Tyreus, B. D.; Luyben, M. L. Plantwide Process Control; McGraw- Hill: New York, 1998. (19) Rase, H. F. Fixed-Bed Reactor Design and Diagnostics; Butterworth Publishers: Woburn, MA, 1990. (20) More, J. J.; Hillstorm, K. E.; Garbow, B. S. MIPACK Project, Argonne National Laboratory, Argonne, IL, Mar 1980. (21) Morbidelli, M.; Varma, A. A generalized criterion for parametric sensitivity application to thermal explosion theory. Chem. Eng. Sci. 1988, 43, 91. (22) Riggs, J. B. Chemical Process Control, 2nd ed.; Ferret Publishing: Lubbock, TX, 2001.

Received for review November 7, 2002 Revised manuscript received March 18, 2003 Accepted April 30, 2003 IE020874P