Ind. Eng. Chem. Res. 1995,34,3872-3877
3872
PROCESS DESIGN AND CONTROL Dynamics and Stability of Polymerization Process Flow Sheets I. Hyanek, J. Zacca, F. Teymour: and W. Harmon Ray* Department of Chemical Engineering, University of Wisconsin, 1415 Johnson Drive, Madison, Wisconsin 53706
In this work, the stability structure of process flow sheets described by systems of differential and algebraic equations (DAE’s)is analyzed. This general framework allows the study of some of the dynamics and stability problems that arise in polymerization processes. A new numerical driver (ISAAC) which is a part of the dynamic flow-sheeting package POLYRED (polymerization reactor design) is used to perform continuation studies and steady-state stability analysis. The Himont liquid pool loop reactor and the Unipol gas-phase fluidized bed reactor are presented as examples of real systems which can become unstable because of the nature of the flow-sheet topology and operating conditions selected.
Introduction
POLYRED Package
The dynamic behavior and stability character of multiunit chemical processes depend not only on the characteristics of the individual process units but also on the topology of the flow sheet. It is well-known that the feedback of energy or mass through heat integration or recycle streams can destabilize a process even when the individual process units are stable. In this paper, the dynamic flow-sheeting package POLYRED, developed by the University of Wisconsin Polymerization Reaction Engineering Laboratory (UWPREL), is used to analyze some of the stability problems that arise in the polymerization processes. Examples are presented of real processes which can become unstable due to the nature of the flow-sheet topology and operating conditions selected. POLYRED is a dynamic flow-sheeting package which provides for the simulation, parameter estimation, optimization, and stability analysis of process flow sheets. The current process unit modules in POLYRED are primarily polymerization reactors, separators, mixers, etc. One of the newest features of POLYRED is the ability to do continuation with respect to module or flow-sheet parameters and to determine the stability character of the process. This allows the precise prediction of operating and design parameters which lead to multiple steady states, oscillatory instabilities, or even chaos. In this paper, these new analysis capabilities will be described and illustrated with real polymerization process flow sheets. Examples will include fluidized bed and loop reactor processes for the production of polyolefins.
The POLYRED Package The POLYRED package is a dynamic flow-sheeting package with the structure shown in Figure 1. There is a CAD environment consisting of a shell around three types of programs: (1)utility routines that allow the user to interactively configure flow sheets, to specify
* Author to whom correspondence
should be addressed.
Present address: Department of Chemical Engineering, Illinois Institute of Technology, Chicago, IL 60616. +
0888-5885/95/2634-3872$09.00/0
CAD ENVIRONMENT
SHELL ,I
I
-.. ....--
Iltllitian
I.,
,I
11
I
General Process Models Kinetic Schemes
1
Drivers Dynamic Simulation
I
Parameter Estimation
Reactor Models Particle Models
u Optimlzation
Figure 1. Schematic structure of the POLYRED package.
parameters and operating conditions for each unit in the flow sheet, to plot the results, etc.; (2) general process models which are the modular units in the flow sheet; (3) numerical driver routines which operate on the flow-sheet model to carry out dynamic or steady-state simulation, parameter estimation, continuation and stability analysis, optimization, etc. This modular structure allows the different programs to be coded and maintained independently so that, for example, new models do not require changes in drivers or utilities and new drivers do not require changes in models. Once configured, flow sheets and models can be stored as a “process” and easily reanalyzed with new parameters, operating conditions, etc. The process models allow polymerization using a choice of five types of chemical mechanisms: free radical, ionic, group transfer, Ziegler-Natta, and condensation. A wide choice of polymerization reactors is available, corresponding to those used by industry. For example, reactor modules include tank, fluidized bed, tube, gas stirred bed, loop, and distillation tower. These reactors can represent multiphase behavior such as suspension, emulsion, dispersion, and slurry. A recent study (Debling et al., 1994) compared the performance of most of the common polyolefin processes when making product grade transitions. This required the dynamic simulation of multistage flow sheets involving multiple
0 1995 American Chemical Society
Ind. Eng. Chem. Res., Vol. 34,No. 11, 1995 3873 CSTRs, multiple fluidized beds, multiple loops and fluidized beds in the same flow sheets, multiple stirred gadsolid reactors, etc. A second study (Appert et al., 1992) analyzed the performance of distillation reactors combined with a train of CSTRs for the transesterification and polymerization stages of a PET process.
Process Stability Analysis One of the newest features of POLYRED is the ability to carry out continuation with respect t o one or more process parameters in order to map the steady-state and dynamic stability behavior of the process in parameter space. The ISAAC (interactive stability analysis and continuation) driver within POLYRED allows the user to accomplish this task interactively without writing any computer code. One needs only to define the flow sheet, provide the parameters for the process units,and specify which parameters are the continuation variables. There are several levels of stability analysis which can be requested. The ISAAC driver analyses the behavior of the set of differential algebraic equations (DAE’s)that arise from the process model. The driver is based on a modified version of the general continuation package, AUTO (Doedel, 1986). Because of the flow-sheet structure, discontinuities, and DAE formulation of the process models, major modifications of AUTO were required (AUTO is designed for only smooth ODE’S). ISAAC provides tools for steady-state continuation and stability analysis of implicit DAE’s. In addition, it is common for POLYRED models to contain IF-ELSE-ENDIF’blocks, leading to discontinuous models. Thus, AUTO had t o be internally reorganized to provide a stepwise mode in addition to the continuous one to allow for the checking of the discontinuities and automatic restart. In addition to following a branch of steady states, ISAAC is capable of the accurate determination of branching points, limit points, and Hopf bifurcation points. Once these are detected, ISAAC can be used to compute the locus of these points in two parameters. Currently, features are being added for periodic solution continuation and stability analysis that will handle the internal model discontinuities mentioned above. The ISAAC driver was written as a fully interactive numerical tool compatible with other POLYRED drivers (simulation, parameter estimation, etc.). Prior to steadystate continuation, the user is allowed to modify process flow-sheet parameters as well as the numerical parameters. The run time screen displays the progress of the continuation together with the information on the special points that were crossed along the branch. All the special points are stored and maintained in the internal ISAAC database. After the continuation is complete, the user is able to browse through the database and display more information on selected points. The available options for the last specified special point are displayed on the screen, allowing the user to select one and resume continuation. The user can save the data to plot files and plot the results in a variety of formats.
Mathematical Formulation The ISAAC driver is based on a modified version of the general continuation package, AUTO (Doedel, 1986). The driver can do a limited bifurcation analysis of process flow-sheet models that are mathematically represented by a system of algebraic-differential equa-
tions of the form
f(x,x’,p)= 0 ,
x, x‘,f E R”,p E R
(1)
Here f a r e residuals of model equations, x and x’ are vectors of state variables and their time derivatives, and p is some flow-sheet parameter. In the first step of bifurcation analysis, the stationary solution branches are determined as solutions (x(s),p(s)) of the nonlinear system of equations f(x,x’=O,p) = 0
(2)
where s denotes some homotopy parameter. Assuming that a solution xo = x(s0) is known for some value po = SO), then the direction (&,bo) of the branch a t this point coincides with the null vector of the n by n 1 dimensional matrix (f,(xo,po>lf,(xo,po)). From the starting point, the stationary solution branch is traced out in a stepwise manner using the pseudoarc-length continuation technique. Assuming that (Xj-1,pj-1) and (kJ-1,,bJ-1) are known, the next solution (xj,pj) can be determined from the equations
+
f(x,x’,p)= 0 AS = eU2(xj -
x~Js,-,
+ eJpj - pj-l)fij-l
(3)
where Ou and 0, are scaling constants corresponding to x and p.
The determination of the local stability of the system of equations (l),as we vary the continuation parameter p, requires the solution of an eigenvalue problem a t the point p = po. To begin the discussion,let us linearize eq 1for a fixed value of the parameter (p = pol
f,dx
+ fd6X’= 0
(4)
where f, and f,? are Jacobians a t the steady state of interest, 6x is a small perturbation about the steady state, and 6x‘ = d(dx)/dt. Then (4) can be written more explicitly as
d(dx) fx.-- -$dx dt Note that, since only some of the components of x are differential states, fx, is of reduced rank. Equation 5 leads directly to the generalized eigenvalue problem
fxv= -Af,,v
(6)
Equation 6 is expressed in terms of the system Jacobian matrixes (fx,fx,),and it does not require any splitting of the system state variables for its calculation. Well-known generalized eigenvalue subroutines (Smith et al., 1976) can then be easily implemented in order t o determined the system stability at any given steady state. The limit points on the steady-state solution branch are determined as the points where the scalar function q(s) = k s ) vanishes. Once a limit point has been detected, the whole branch of limit points in two parameters p and 17 can be obtained by applying the pseudoarc-length continuation technique to the enlarged system of equations comprised of the system steadystate equations, the condition of a zero eigenvalue applied to eq 6 , and a normalizing eigenvector condition
3874 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 Purge
Purge
Catalyat Feeder
800 Te 700 (K) 600 Ethylene Comonomer Hydrogen
Ethylane Comonomer Hydrogen
A,kvl
1
No Catalyst Decay 1 - Low Hydrogen Level 2 - High Hydrogen Level ..
..... .......... Me1 ting ....... ................ ::............................... .......
nlkvl
....__...
,_.'
Product
Product
Figure 2. Fluidized bed reactor process for the production of polyolefins.
(required for unique eigenvector calculation)
f(x,x'=O,p,v)= 0
(7)
fx(x,x'=0,p,1;I)v = 0
(8)
vTv - 1 = 0
(9)
+
The previous system is composed of 2n 1equations in 2n 2 variables, (x, v, p, a). The detection of the Hopf bifurcation points is accomplished by monitoring of the sign of the scalar function q(s) = Re{I(s)}, where Im{L(s)} = w > 0 and I is the complex eigenvalue with minimum imaginary part. This implies that, a t a Hopf bifurcation point, there is a purely imaginary eigenvalue of the form I = io to which corresponds an imaginary eigenvector of the form v = Re{v} iIm{v}. As in the case of turning points, in order t o uniquely determine the eigenvector, a normalizing condition is required. Applying these conditions to eq 6, the loci of Hopf bifurcation points in two parameters are obtained by continuation of the solutions of the enlarged system
+
+
(Re{v})TRe{v}
+ (Im{v})TIm{v}= 1
(11)
(Im{v})TRe{v} - (Re{v})TIm{v} = 0 (12) This is a system of 3n 2 equations in 3n 3 variables, (x,Re{v},Im{v}, w , p, 7).
+
+
Some Examples Although any flow sheet can be analyzed by ISAAC, space limitations require us to limit our discussion to two types of processes here. First, we will analyze the behavior of single-stage and multiple-stage fluidized bed processes for the production of polyolefins. The flow sheet for the two-stage process is shown in Figure 2. For the specific example to be considered here, we will be making a linear low-density polyethylene (LLDPE), which is a copolymer of ethylene with an a-olefin such as propylene, butene, hexene, etc. Hydrogen is added to the reactor as a chain-transfer agent t o control the molecular weight. LLDPE can be made in a single reactor, or for a bimodal molecular weight distribution (MWD), the product will be made in two stages as shown in Figure 2. One has the option of making the high
,..__..... -...*...."
i
Ind. Eng. Chem. Res., Vol. 34,No. 11, 1995 3875 1000
1
n
1000
1
r
Te ' 0 °
Te 700
(K) 600
400 300
High Hydrogen Level
1
(K) 600
500
1 1
Melting
.....................................................
400
300 0.00
...................
0.02
0 04
0.06
0.08
H Figure 7. Steady-state continuation showing stage 2 for two-stage 0
4000 6000 8000 10000 Time (sec) Figure 4. Reactor runaway trajectory after temperature controller failure (without catalyst decay): T f = 300 K T, = 343 K, u/umt = 10.0; high H2 level.
2000
process without catalyst decay: (0)Hopf bifurcation; ( 0 )operating steady state (T, = 368 K u/u,f = 10.0; Tf = 320 K).
520 1
(K) 600
Te 8oo
'$00
I
I !
............................................................................
200
320 300
Time (sec)
Catalyst Decay 1 - Low Hydrogen Level 2 - High Hydrogen Level
340
' 0
Figure 8. Reactor runaway trajectory after simultaneous temperature controller failture (without catalyst decay). I
0.4 0.6 0.8 1 qcat (41s) Figure 5. Steady-state continuation with catalyst decay: 1,low H2; 2, high Hz; ( 0 )Hopf bifurcation; ( 0 )operating steady state (T, = 333 K u/umf = 7.0; Tf= 300 K).
480
0.2
460 440
r
I
i Melting
Te 400 (K) 380
1 - First Reactor 2 - Second Reactor
500
1
480
/
1
.....
2000 4000 6000 8000 10000 12000
0
1
460
320 300 0
440 Te
0.10
qcat (gk)
420
0.2
0.4 0.6 qcat (ds)
0.8
1
Figure 9. Steady-state continuation for a two-stage process with catalyst decay: (0) Hopf bifurcation; ( 0 )operating steady state.
400 (K) 380 360
340 320 300 0
20000
40000
60000
80000 100000
Time (sec) Figure 6. Reactor runaway trajectory after temperature controller failure (with catalyst decay): Tf = 300 K, T, = 333 K, u/umf= 7.0; high H2 level.
catalyst feed rate in the case of a catalyst having no significant decay. Because of the inhibiting effect of H2 on the catalyst activity, much more catalyst is required for operation at high Hz levels. Note that the steadystate operating temperature of 343 K requires operation beyond the Hopf bifurcation point in both reactors so that the steady state is open loop unstable. With sufficient cooling capacity, this steady state is easily stabilized by a PID temperature controller. Unfortunately, if this controller fails, the reactor will run away in a few hours, as shown in Figure 4. Without countermeasures (catalyst poisoning, etc.), this can result in meltdown of the polymer powder into a large molten mass of polyolefin inside the reactor.
In the case of a catalyst with significant catalyst decay, the bifurcation structure is shown in Figure 5. Note that the accelerating catalyst decay with increasing reactor temperature has eliminated the multiple steady states shown in Figure 3. However, the unique branch of steady states still has Hopf bifurcation points so that temperatures in the normal steady-state operating range (333-373 K) are open-loop-unstableoperating points. These are easily stabilized by sufficient cooling capacity and a PID temperature controller. In this case, controller failure will lead to sustained oscillatory behavior, as shown in Figure 6. These oscillations could be survived by reactor except that the temperature peaks exceed the melting point of the polymer by more than 50 K. Thus, the reactor will form some molten polymer which will create large chunks of polymer in the fluidized bed, possibly finally leading to meltdown. The specific details of the stability behavior and the resulting dynamics depend strongly on the type of catalyst and the reactor operating conditions. The use of the ISAAC driver offers the possibility that, by suitable catalyst formulation and choice of steady-state operating conditions, one could design the process to be much more robust to controller failure and potential run
3876 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 550
I
I
500
Te
450
(K) 400
350
___
300
0
50000
100000
Time (sec)
150000
200000
Figure 10. Reactor 1 dynamics after simultaneous temperature controller failure (with catalyst decay). 420 I
0
I
50000
100000
150000
I 200000
Time (sec) Figure 11. Reactor 2 dynamics after simultaneous temperature controller failure (with catalyst decay).
away-for example, by making catalyst decay more temperature sensitive so that the oscillations in Figure 6 remain below the polymer melting point.
Multiple-Stage FBR Operation The basic operating conditions for the multistage process for the production of a bimodal MWD are also given in Table 1. Here we will examine the stability of the two-stage process when we wish to produce a bimodal MWD by making low molecular weight polymer in stage 1(high H2) and high molecular weight polymer in stage 2 (low H2). As in the case of single-stage reactors, we carry out the analysis for a catalyst without deactivation and then for a catalyst having significant deactivation. Figure 7 shows the steady-state continuation diagram for the second stage of the process for a catalyst without
decay. The first-stage behavior is shown as curve 2 in Figure 3. Note that, in stage 2, the steady-state behavior is much more complex than stage 1. There is even a region of nine steady states, although most of them are above the polymer melting point. The simulations show that stage 1 is stable while stage 2 is unstable for the catalyst feed rate of qcat = 0.03 gls. For simultaneous temperature controller failure in the process, Figure 8 shows that reactor 2 runs away while reactor 1remains stable, as expected. Thus, stable twostage operation requires a balancing of reactor sizes, diluent concentration, and catalyst feed rate. When there is catalyst deactivation, the two-stage process operates much better because the low H2 levels in stage 2 are balanced by decreased catalyst activity due to deactivation. The continuation diagrams in Figure 9 show that both stages would operate at an unstable steady state so that simultaneous temperature controller failure leads to oscillations in each reactor (cf. Figures 10 and 11). Note that the oscillations in stage 2 are more complex initially because the naturally occurring autonomous oscillations are also driven by oscillations in the feed conditions. However, after some time, the stage 2 oscillations synchronize to the same period as reactor 1. It is interesting that the stage 2 oscillations remain below the polymer melting point so that the reactor could operate if the stage 1 polymer coming as feed did not have too many large chunks.
Single-Loop Reactor As a second example process, we consider the behavior of loop reactors for the polymerization of polyolefins. Let us analyze the behavior of the lead loop reactor in the process for the polymerization of propylene shown in Figure 12. This technology was developed by Himont to produce polypropylene homopolymer and impact polypropylene. The details of the process and the loop reactor model may be found in Debling et al. (1994) and Zacca and Ray (1993). For the process conditions and catalyst system described in Zacca and Ray (1993), the steady-state continuation diagram of loop exit temperature with respect to loop recycle ratio is shown in Figure 13. Spontaneous catalyst decay is assumed using representative parameters (e.g., Hutchinson (1990) and Drusco and Rinaldi (1984)). Note that the reactor is unstable over a wide range of recycle ratio values,
I
Purge
t &
Hydrown
cw Prowlene ._ I
CocatalYat ( AI - Alkyl )
I
m r
N2
Prepolymerlzatlon
Finishing Loop Reacton 7OoC , 3 6 40 atm
-
Figure 12. Multistage loop reactor process for the polymerization of propylene.
Fluidized Bed Reactor EOOC, 20 atm
Ind. Eng. Chem. Res., Vol. 34,No. 11, 1995 3877 175
To
of the University of Wisconsin Polymerization Reaction 7 Engineering Laboratory (UWPREL) for support of this
:\
150 125
(deg C)
research.
Nomenclature Critical Temperature
100
‘
50 2.5
I
5
7.5
10
12.5, 15
17.5 20
Recycle Ratio
Figure 13. Steady-state continuation for the loop reactor; loop exit temperature vs the recycle ratio: (0) Hopf bifurcation; ( 0 ) operating steady state. 200 7
“V
0
50000
100000
150000
Time (sec)
Figure 14. Reactor exit temperature oscillations for a recycle ratio of 15.
and the small region of multiple steady states is completely within the range of unstable steady states. Thus, the reactor oscillates for all values of recycle ratios within this range. As an example, for a recycle ratio of 15, the reactor exit temperature oscillates as shown in Figure 14. Note that the exit temperature peak exceeds the melting point of polypropylene and also rises above the critical temperature of propylene by more than 100 “C so that there could be reactor fouling from molten polymer and “vaporization” of propylene. This could lead to traveling waves, mechanical stress, and process shutdown. The actual consequences would depend on a more detailed analysis. The composition and temperature profiles in the loop for this and other flow-sheet configurations will be presented in a future publication.
Conclusions The use of POLYRED for stability analysis of polymerization process flow sheets described by a system of differential and algebraic modeling equations (DAE’s) has been illustrated with a number of examples. This feature, unique in a process simulation package, is shown to be a very useful tool in developing design, operating, and control procedures to ensure productive and safe process operation.
Acknowledgment We are indebted to the National Science Foundation, CNP,-Conselho Nacional de Desenvolvimento Cientifico e TecnoMgico-Brasil, and the industrial sponsors
f = vector of residuals of the modeling equations fx = matrix of derivatives off with respect to x fx,= matrix of derivatives off with respect to x’ qcat= catalyst mass flow rate (g/s) QF = feed volumetric flow rate (cm3/s) QR = recycle volumetric flow rate (cm3/s) Rec = QR/QF,recycle ratio (dimensionless) s = homotopy parameter Tf= feed temperature (K) T, = operating (steady-state) temperature (K) u = inlet gas velocity ( c d s ) umf = minimum fluidization velocity ( c d s ) v = eigenvector Re{v}, Im{v} = eigenvector real and imaginary parts x = state vector x’ = time derivatives of the state vector S = state vector derivatives with respect to the homotopy parameter Greek Symbols il = vector of eigenvalues p , 7 = flow-sheet parameters p = flow-sheet parameter derivative with respect to the
homotopy parameter
Literature Cited Appert, T.; Jacobsen, L.; Ray, W. H. A General Framework for Modeling Polycondensation Processes. Proc. 4th Znt. Workshop React. Eng., DECHEMA Momgr. 1992,127, 189-197. Choi, K.; Ray, W. H. The Dynamic Behaviour of Fluidized Bed Reactors for Solid Catalysed Gas Phase Olefin Polymerization. Chem. Eng. Sei. 1986,40,2261-2279. Debling, J. A.; Han, G. C.; Kuijpers, F.; Verburg, J.; Zacca, J. J.; Ray, W. H. Dynamic Modeling of Product Grade Transitions for Olefin Polymerization Processes. MChE J . 1994, 40 (31, 506-520. Doedel, E, AUTO: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations. Technical Report, Caltech, 1986. Drusco, G.; Rinaldi, R. Polypropylene-Process Selection Criteria. Hydrocarbon Process. 1984, 113-117. Hutchinson, R. Modelling of Particle Growth in Heterogeneous Catalyzed Olefin Polymerization. Ph.D. Thesis, University of Wisconsin-Madison, 1990. Smith, B. T.; Boyle, J. M.; Dongarra, J. J.; Garbow, B. S.; Ikebe, Y.; Klema, V. C.; Moler, C. B. Matrix Eigensystem RoutinesEZSPACK Guide; Springer: Berlin, 1976. Zacca, J. J.; Ray, W. H. Modelling of the Liquid Phase Polymerization of Olefins in Loop Reactors. Chem. Eng. Sei. 1993,48 (221, 3743-3765.
Received for review March 13, 1995 Revised manuscript received J u n e 16, 1995 Accepted J u p e 23, 1995* IE950176N
Abstract published in Advance ACS Abstracts, September 15, 1995.
