244
Ind. Eng. Chem. Res. 1992,31, 244-254
Design and Analysis of Chemical Processes through DYNSIM Rafiqul Gani,* Esben L. Smensen, and Jens Perregaard Engineering Research Centre IVC-SEP, Chemical Engineering Department, The Technical University of Denmark, 2800 Lyngby, Denmark
This paper shows how an equation-oriented process simulator can be employed to solve simulation problems related to the design and analysis of chemical processes. The flexibility of equation-oriented process simulators allow the development and use of problem-specific simulation strategies by considering the particular needs of the specified design/analysis problem. These strategies allow switching between steady-state and dynamic modes of simulation, switching between stiff and nonstiff methods of integration, and switching between modes of solution (solving a set of ordinary differential equations plus procedures or a set of differential-algebraic equations) for the specified problem. The integration of simulation and design/analysis of chemical processes is illustrated through four practical examples.
Introduction Dynamic and steady-state simulationsplay an important role in the design and analysis of chemical processes. Through dynamic simulations, the effect of changes in design variables on the continuous operation of a plant can be analyzed. Dynamic simulations can also be employed for studying the effect of disturbances on the transient behavior of the process. Steady-state simulations provide information on the operating conditions of a process and the effect of disturbances on the process stream variable(s). Reliable design and the correct analysis of a chemical process need, therefore, correct prediction of the behavior of the process. Often, chemical processes contain highly nonideal mixtures whose physical properties are sensitive to the design and operation of the process. Accurate simulation of these processes needs to employ appropriate thermodynamic models. Usually, these thermodynamic models are complex and require special handling in order to keep the computational time at acceptable levels. For example, use of thermodynamic models combining activity coefficient models like UNIFAC (Fredenslund et al., 1977) with equations of state are recommended for simulation of supercritical/near-criticalextraction processes. Equation-oriented process simulators offer the flexibility needed for adapting the simulator for various types of simulation problems related to the design/analysis of chemical processes. For example, they offer easy integration of steady-state and dynamic simulation during design/analysis of chemical processes. Also,since all the features of the simulator may not be needed for a specific simulation problem, the equation-oriented process simulators allow the possibility of employing problem-specific simulation strategies. For some design/analysis-related simulation problems, a combination of dynamic simulation (nonstiff integration method) and steady-state simulation provides the needed information. For another problem, a combination of nonstiff and stiff methods of integration is needed for the dynamic simulation part. Also, the dynamic simulation can be performed in the DAEs (solving a set of differential-algebraic equations) mode or in the ODES (solving a set of ordinary differential equations plus procedures for algebraic equations) mode. For some simulation problems, the nonlinear algebraic equations resulting from the prediction of phase equilibria can be linearized (or simplified), while, for other problems, linearization may lead to significant errors. Recently, Gani et al. (1990) have shown how problemspecific simulation strategies can be developed in an equation-oriented process simulator. Gani et al. (1990) used the equation-oriented dynamic simulator, DYNSIM
(Sarensen et al., 1990),as an example. DYNSIM is a mixedmode dynamic process simulator which can also solve some steady-state simulation problems for a broad range of chemical proteases. It is mixed mode because it allow flow sheets with dynamic as well as nondynamic units, dynamic and steady-state modes of simulation, and ODE and DAE modes of solution. One of the aims for developing DYNSIM is that it should be particularly suitable for problems related to the design and analysis of parts of chemical processes containing highly nonideal mixtures: for example, the separation unit or the reactor/separator unit of a chemical process containing nonideal mixtures or operations at extreme conditions of temperature and pressure. Gani et al. (1990) have demonstrated some of the predictive capabilities of DYNSIM for different types of simulation problems. The objective of this paper is to extend the problemspecific simulation strategies of Gani et al. (1990) together with the special modeling/computational featurea that can be employed in an equation-oriented simulator and to illustrate how these features can be employed in integrating process simulation with process design and analysis. This paper therefore describes the problem formulation, the method of solution, and the results obtained through DYNSIM for different simulation problems related to the design/analysis of chemical processes.
Problem Formulation The design and analysis problems being considered in this paper are defined as follows: the study of the behavior of chemical processes due to changes/disturbances in equipment sizing parameters (a,) and/or manipulative variables (d2),determination of the new values of the process stream/unit variables (p), and also, the study of the effect of specified values of p on dl and/or d2. Equipment-sizing parameters are those related to the size of the equipment. For example, in a flash unit they are cross-sectional area, volume, and valve coefficients. The manipulative variables are those that are usually specified in open-loop (without control) simulations or manipulated during closed-loop simulations. Typically, they include variables like reflux rate, vapor boil-up rate, reboiler heat duty, etc., in a distillation column. The process stream/unit variables are those related to the output streams of a process/unit. Variables like percentage recovery, product yield, flash tank pressure, stream enthalpy, etc. are examples of process stream/unit variables. For dynamic simulation,all variables belonging to dl and d2are specified and all elements of p are calculated. The size of the arrays of dl and d2 also reflect the number of
O888-5885/92/2631-O244$03.OO/O 0 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 245 degrees of freedom. Some steady-state simulation problems specify all variables from dl and dF Other steadystate simulation problems may specify variables from p with the condition that if n variables from p are specified, k-n variables from d, and/or d2are left unspecified and are therefore calculated. k is the dimension of the vector d consisting of all elements of dl and dP The simulation problems related to the design and analysis of chemical processes are therefore divided into the following three types (given, a reference operating condition or a default initial condition): (i) the study of the effect of changes in dl on the transient behavior (dynamic simulation) and the new operating condition (steady-state simulation); (ii) the study of the effect of disturbances in d2 on the transient behavior (dynamic simulation) and the new operating condition (steady-state simulation); (iii) the study of the effect of changes in p on the unspecified d, and/or dz (steady-state design problem).
equipment-sizing parameters on the conditions of operation of the process with both modes of simulation. The set of implicit algebraic variables, x, usually representa the variables associated with computations of phase equilibria. For simulation problems of type iii, some elements of p are specified instead of the same number of variables represented by d. For example, specification of tank liquid level height (belongs to p2) instead of the tank cross-sectional area (belongs to d,) provides the value of the area needed to match a specified liquid level. Typically, eq 2 represents the subset of algebraic equations which may be solved separately as procedures or added to the ODES to be solved in the DAE solution mode. An example is the equations subset representing the prediction of phase equilibria.
+ K i m i- n:)nF/Cni v = C(uF + uy,
(5)
Simulation Model The dynamic and steady-state simulator, DYNSIM (Scarensen et al., 1990),has been used in this paper for all simulations related to design and analysis of chemical processes. DYNSIM has been chosen because it has been developed at the Institute for Kemiteknik, The Technical University of Denmark, Lyngby, Denmark, and because it has a highly structured and flexible simulation strategy which allows one to take advantage of the particular needs of a specific simulation problem. In principle, any other equation-oriented simulator could have been used. A brief description of DYNSIM, highlighting the special features developed and/or used in this work, is given below. Modeling Aspects. The models for all the equipment presently available in DYNSIM are derived such that the equations form a differential-algebraic equations (DAEs) system for dynamic simulations and correspond to an INDEX1 problem (as defined by Pantelides et al., 1988) or form an algebraic system (AE) for steady-state simulations. It should be noted that the same set of equations is used for either mode of simulation. It is beyond the scope of this paper to describe the details of the different unit operation models available in DYNSIM. These are given by Scarensen (1991) and will be made available to the reader upon request. Very briefly, the model equations consist of first-order ordinary differential equations (ODES)representing the mass and energy balance equations (eq l),a set of implicit algebraic equations (IAEs, eq. 2), and, a set of explicit algebraic equations ( E m s , eqs 3 and 4). By setting D D = dy/dt = f(y,x,p,d,t) (1)
T = f(hL,P,nL) u t = f(T,P,nL)
(7)
uy = f(T,P,nV)
(9)
x = f(y,x,d) (2) PI = f(Y,d) (3) ~2 = fb, ,y,d) (4) = 0 in eq 1,the set of AEs for steady-state simulation is obtained. The IAEshave more than one unknown variable while the E m s have only one unknown variable. The vector d represents all elements of dl and d2 while the vector p represents all elements of p1 and pz. For all simulation problems of types i and ii, the vector d is specified along with the initial condition for y and the initial estimates for x. In this work, for dynamic simulations, y represents the set of differential variables, while, for steady-state simulation, y represents the set of unknown algebraic variables. In both modes of simulation however, y represents the mass and energy accumulation terms. This makes it possible to study the effect of
ni = nt
(6) (8)
Ki = f(TP,n,nL)
(10) Equations 5-7 represent the IAEs (eq 2), while eqs 8-10 represent the EAES (eq 4). The implicit algebraic variables are nL,T,and P, while the explicit algebraic variables are vL,v", and K. n is a differential variable belonging to the vector y and V is a specified variable (in this case, the tank volume) belonging to the vector d. Typically eq 3 relates the flow parameters of a stream to the differential and specified variables. For example, the flow rate of a stream coming out of a tank is computed from
F = f(n,Rm)
(11) In eq 11,n are the known differential variables while R" is the specified design variable (valve coefficient). A similar equation for a flash tank, however, belongs to eq 3.
FL= f(T,P,RL,nL) (12) In eq 12, T,P, and nL are implicit algebraic variables belonging to the vector x. Another example of equations of type 4 is the equations calculating some property (for example, enthalpy, density, etc.).
hV = f(n,nL,T,P) (13) In the above equation, h", the vapor enthalpy, is only calculated after nL,T,and P have been determined. Thus, while eqs 8-10 are repeated for every function evaluation for the IAEs, eq 13 is solved only once per iteration. As can be seen from above, the nature of eqs 2-4 requires a predefined calculation order for the evaluation or solution of the right-hand sides (RHSs) or solution of these equations. This classification and ordering of equations make it easy to integrate process simulation and design through problem-specific simulation strategies. Simulation Strategy. The grouping of the model equations into the three types, as shown above, makes it possible to develop a highly flexible and structured simulation strategy, which can be adapted for different simulation problems. The main features of this simulation strategy is described below. ODE Mode vs DAE Mode. If only eq 1 is solved simultaneously by the integration method and eqs 2-4 are solved separately, as procedures, we have the ODE simulation mode. If eq 2 is included in the set of equations to
246 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
be solved simultaneously, we have the DAE simulation mode. Therefore, for the ODE mode, the calculation order is eq 3, eq 2, eq 4, and finally the RHS of eq 1, while for the DAE mode, the calculation order is eq 4, eq 3, and the RHS of eqs 1and 2 (note that in this case values for both x and y are known at any time t or iteration k). Equations 3 and 4 are not included in the DAE set of equations since they are explicit equations with only one unknown variable. Also, they do not need any iterations for their solution. Some equations (for example, eq 10) may in some cases (depends on the thermodynamic model) need iterations. From a computational point of view and our experience, it is better not to include them in the DAE set. For both modes of solution, a consistent initialization of x and y is needed. DYNSIM provides this by first determining y, x, and unspecified elements of p and d from the specified values of p and d (eq 3 for y, eq 2 for x, and eq 4 for unspecified p or d for a type iii problem) and then starting all dynamic simulations with the ODE mode. After a specified number of integration steps, a switch to the DAE mode is made. Note that this procedure is only needed if the reference steady-state condition or the initial condition is not known. Since the procedures for solving the IAEs in the ODE mode are more robust than the method of solution in the DAE modes, it is possible to start simulations from a default initial condition which does not satisfy eq 2 (consistent initialization in this paper includes a balance between satisfying eqs 1-4 and determining physically meaningful values for the stream/unit variables). For transient responses due to disturbances around a reference steady state (problems of types i and ii), switching after the first integration step is recommended. For larger changes or for simulations starting from a condition far from the steady state, the transient behavior needs to be considered when the time for switching is being determined. Generally, the computational time for the DAE mode is smaller than the ODE mode. By solving different simulation problems, we have noted that the computational times with the two modes of simulation vary with the type of units in the process flow sheet. For example, for a flash unit, the DAE mode needs only about 20% of the computational time needed by the ODE mode (Gani et al., 1990), while, for distillation columns, the computational time varies with the type of simulation problem (Perregaard et al., 1991). The reason for this variation can be attributed to the efficiency of the method of solution employed to solve the algebraic equations separately as procedures, the thermodynamic model used, and the number of additional equations in the DAE mode. Also, for some unit operation models, there is no significant difference between the ODE and DAE modes. This is because the number of equations being added are small and they are linear or very close to being h e a r (for example, tank mixer, simple continuous stirred tank reactor (CSTR) reactor, etc.). For a flash unit (closed tank), to solve the EAEs and IAEs as procedures, a three-loop iterative method is needed. That is, for a given total molar holdup (represented by y) and the flash unit sizing parameters (a),the temperature (T), the pressure (P),the equilibrium-phase compositions, and the ratio of vapor/liquid molar holdup are determined. Thus, since the number of additional equations is not very large, even though their percentage increase is large, computational time with the DAE mode is saved by avoiding the three-loop-iteration method. It should be noted that the procedures employed by DYNSIM are very robust (nearly never fails to converge, even from
bad initial estimates) and reliable. For diatillation calculations however, the set of eqs 2 and 3 is much easier to solve as procedures (needing only a one-loop-iteration method) and the time saved in the DAE mode by avoiding the iterative method is counterbalanced by the increase in the number of equations. This is because the procedures in the ODE mode achieve convergence very rapidly (maximum of one or two iterations are needed during the numerical evaluation of the Jacobian matrix). The size of the Jacobian matrix in the DAE mode however is nearly twice as large compared to the ODE mode. Dense, Sparse, and Partitioned Modes. For many chemical processes, the Jacobian matrix J = dfi/dyj i = 1(1)N,j = l(1)N (14) needed for the “stiff“ integration method and for the Newton-like method employed for steady-state simulation, can have many zero elements. The advantage of this sparsity of the Jacobian matrix is taken through the use of sparse techniques (Duff, 1979) or through a technique for partitioning the equation set into subsets and solving all the partitions simultaneously (Sarensen, 1991). The advantage of the “partitioned” mode (when applicable) is that different integration methods can be employed for different partitions and, as in the sparse mode, the number of function evaluations needed to evaluate the Jacobian matrices are drastically reduced. Note that this “partitioned” mode is different from the independent/ parallel modular techniques proposed by Hillestad and Hertzberg (1986) because it does not need interpolations/exhapolations or iterations on a time horizon for flow sheets with recycle streams. Dynamic vs Steady-State Mode. The same set of equations (eqs 1-4) are solved for steady-state and dynamic simulation modes. For problems of types i and ii, the set of differential variables (y) and state variables (x)are also the same. Thus, for these problems, a switch between the two modes of simulation can be made after any integration step (or iteration). The switch from the dynamic mode to the steady-state mode is however recommended only for simulationsaround a reference steady state with respect to problems related to determining the reference steady state (when it is unknown) or to studying the effect of changes/disturbances on a process operating at the reference steady state. This limitation is due to the need for good initial estimates for Newton-like methods. Further details of this switching policy can be found in Gani and Cameron (1989). For problems of type iii, since some elements of p are specified instead of d, the corresponding EAEs are added to the set of IAEs and solved either as procedures or as a set of AEs. Computational Aspects. The structure of DYNSIM is shown in Figure 1. DYNSET and DYNDAT are two utility programs related to DYNSIM. The program DYNSET determines the size of the different arrays from the flow-sheet information. The program DYNDAT generates the input files for DYNSIM. DYNDAT contains a knowledge base of groups, model parameter tables, parameters for the numerical method, and information on the units and their associated design/manipulative variables. Presently, DYNSIM has the following thermodynamic models for computation of phase equilibria: nine equations of state, five activity coefficient models, the GCEOS model (SkjoldJargensen, 1988) and MHV2 (Dah1 et al., 1991). Also, DYNSIM has process unit models for diatillation columns, flash tanks, CSTR reactors, plug flow reactors, heat exchangers, stream divider, evaporators, condensers, decanters, and different types of mixers. The numerical routines
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 247 Table I. Simulation Strategies for Different Problems' simulation mode at various times ( t )
>>>>
>>>>I
TERMINAL DYNSIK.RES DYNSIK.KIP
1
problem type i, ii iii
to
tl
t2
t3
1,1,1,2,2 1,1,1,1,2 2,1,1,1,2 2,1,2,1,2 2,1,1,1,1 2,2,1,1,1 2,2,2,1,1 2,2,1,1,1 2.2.2,l.l simulation mode at various times (t)tr
tn
tn
'The five parameters have the following meaning: (parameter 1) mass balance = 1, energy balance = 2; (parameter 2) ODE mode = 1, DAE mode = 2; (parameter 3) dynamic mode = 1, steadystate mode = 2; (parameter 4)stiff method = 1, nonstiff method = 2; (parameter 5) model simplification/reduction = 1, full model = 2. Two seta of values at the same time indicate two different possibilities. Figure 1. Structure of
DYNSIM.
(Cameron and Gani, 1988) include "stiff" (backward difference formula (BDF) method and diagonally implicit Runge-Kutta (DIRK) method) and "nonstiff" (AdamsMoulton and Runge-Kutta) methods for integration. The purpose of subroutine ADMIN (see Figure 1)is to communicate between the integrator and the process unit models and to determine the problem specific simulation strategies. DIFFUN evaluates the rightrhand side of eq 1(and that of eq 2 in the DAE mode). For all process units, only mass balance or mass and energy balance can be performed. Problem-SpecificSimulation Strategies for Design and Analysis. Because of the nature of simulation problems related to design and analysis, different aspects of the simulation strategy described above may not be needed for a given simulation problem. For simulation problems i and ii, we propose converting some of the IAES into EAEs through model simplification. For example, for distillation columns, the IAEs represent the equations related to the computation of phase equilibria. The phase i(l)NC, which is given by eq 10 equilibrium constant, Ki, can be simplified to the following form by regression with data from a reference steady-state condition: In Ki = ai b i / T (15)
+
In the above equations, NC is the total number of compounds present in the system and subscript i indicates compound i. The total number of data points for each fit is equal to the total number of plates in the distillation column, thereby providing a range of application with respect to temperature and mixture composition. The coefficients ai and bi are applicable only for simulations around the reference condition which is used for fitting aiand bi. For simulation problems where only the exit conditions are needed, model reduction rather than model simplification is proposed (Gani et al., 1990). If separation process unite are present in the process flow sheet, then they are ideal candidates for model simplification/reduction. In this paper, under problems of type iii we include design problems where a stream variable(s) is specified instead of a unit parameter(s). If the specified stream variable (belonging to p) and the unknown unit parameter (belonging to d) belong to the same unit, then the proposed simulation strategy simply rearranges the equations and unknown variables related to the EAEs since, for these problems, only the EAES are affected. For problems where the specified stream variable(s) and the unknown unit parameter(s) do not belong to the same unit, an extra set of equations are added. These additional equations are
usually related to the equations for a controller (for example, a PI controller). It should be noted that arbitrary exchange of specified variables between p and d are not allowed because it may lead to a singular set of IAEs and/or EAEs. The simulation strategies for problems of types i-iii are given in Table I. The simulation strategies are described through the values of five parameters at different intervals of time (to,tl, ..., t6). For example, for problems of types i and ii, if the reference steady state is not known, the simulation is started from a default initial condition with only the mass balance option, ODE mode, dynamic simulation with the nonstiff option, and use of the full model. Since, initially, the reference steady state is not known, from our experience, it is better to obtain a "good" composition profile with only mass balance before including the effecta of temperatures and pressures. Also, since the default initial estimates for x may not satisfy eq 2, it is safer to start with the ODE mode. The reason for using the nonstiff option for dynamic simulation is that, initially, the step size is small enough to make the nonstiff method stable. At time tl, a switch to the stiff option is made because now, for the step size, the nonstiff method is no longer stable. At time t2,an energy balance is added because a good profile for the composition has been obtained. At time t,, a switch to the steady-state mode is made because the initial estimates are acceptable for the Newton-like method for the steady-state solver and the reference steady state is obtained. From time t4,studies of the effect of changes/disturbances around the reference steady state are made by changing the value of a variable, integrating (dynamic simulation) until time t6 (that is until part of the transient behavior needed to evaluate the time constants have been obtained), and then switching to the steady-state option to determine the new steady state (that is, determine the steady-state process gains). At time t4, a switch to the simplified model is also made if the changes are within * 5 % . Also, at time t4,a switch from the ODE mode to the DAE mode is made if model analysis suggests that computational time can be saved. For problems where the reference steady state is known, it is also possible to start directly from t4instead of to. For problems of type iii, two problems are solved. In the first problem, for a given set of values of p and d, the unspecified p and d together with the unknown x and y are determined (this problem is solved to determine the default initial conditions and to determine the values of y and/or x for a specified p). In the second problem, it is assumed that the reference steady state is known and the effect of changes in design is studied by specifying some element of p and
248 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 Table 11. Designed Unit Parameters for Problems 1 and 2
Feed Stream mole fraction x i feed no.
F /(kmol/ h)
T/K
1
200
1
100
x2
XR
Xd
XK
330
P/atm Xl Cavett Process 1.5 0.05
0.18
0.27
0.30
0.20
270
l-Hexene Process 10.0 0.50
0.00
0.50
Mixer Units mixer no.
V/m3
A/m2
R"
0.20 0.20 0.10
400.0 440.0 330.0
Cavett Process 0.20 0.20 0.10
1 2
3
l-Hexene Process 10.0
1
50.00
4.0 Flash Units
flash no.
A/m2
RL
2 3 4
0.14 0.30 0.18 0.12
0.18 0.80 0.20 0.16
Cavett Process 320 300 260 230
1 2
0.50 0.50
0.50 0.50
l-Hexene Process 0.0 60
1
variable V/m3 A/m2
V/m3
value 0.50 5.60
R"
Pmh/atm
Q/(MJ/h)
75 35 85 80
0.80 0.80 0.80 0.80
1335.0 -1650.0 3524.0 1020.0
30 30
1.00 1.00
Plug Flow Reactor Unit: Only for the l-Hexene Process variable value variable value variable E/(MJ/kmol) 290*o T&ll/K 0.95 po 1.38 x lo8 HR/(MJ/kmol) U (MJ/m2 K-' h-' ) 100.0 ko
determining the corresponding unknown element of d. The simulation strategy for this problem is given in Table I. Two values in a column indicate two feasible strategies. The simulations start with mass and energy balance (because the reference steady state is known), ODE mode (because the initial estimate may not be good enough for the DAE mode), dynamic simulation or steady-state simulation (depending on what is desirable), the stiff method (because the problem is most probably stiff), and the simplified model (since model simplification is feasible). After the steady-state solution has been obtained at time t2,a switch to the dynamic model is made for a new run with a different problem specification a t time t3. Alternatively, it is possible to obtain the new steady state directly by the steady-state simulation mode.
Application of Dynsim for Design and Analysis Four chemical processes of varying degreea of complexity have been selected for different simulation problems related to design and analysis. The first process is a modification of the well-known Cavett problem (Cavett, 1963). The process flow sheet has extra feed streams, and a five-component (butane, pentane, hexane, heptane, and octane) mixture is considered. In the second process, a simplified flow sheet for the production of l-hexene from l-propene is considered. In the third process, the highly interactive aromatic extraction process is analyzed, and, in the fourth process, the model simplification feature for problems of types i and ii is tested for a near-critical extraction (NCE) process. One interesting aspect of these processes from a simulation and, therefore, design, point of view is that different classes of thermodynamic models are needed for correct the design/analysis of these processes. For example, for the first two processes, use of an equation of state is ap-
value 254.0 82.5
variable T,-/K T,,,/K
value 300.0 303.0
propriate, while, for the aromatic extraction process, a GE-based model is needed for computation of the liquidphase activity coefficients and, for the NCE process, models like GCEOS (Skjold-Jrargensen, 1988) or MHV2 (Dahlet al., 1991) are most appropriate. Note that for the aromatic extraction process, the two-model approach (yep) needs to be employed for the computation of vapor-liquid equilibria (VLE). In this work, the UNIFAC (VLE) model with parameters from Tiegs et al. (1987) is used. The Modified Cavett Process. The Cavett process has four interlinked flash tanks operating at various levels of temperatures and pressures. The specified temperatures and pressures make the flash tanks operate very close to the one-phase region. The flow sheet for the modified Cavett problem is shown in Figure 2a and the specified values for the design variables are given in Table II. Note that, besides these variables, the only other variables needed to start the simulations are the initial condition (or estimate) for x and y. DYNSIM can generate default values for these variables. Thus, the final steady state reached with other simulators should be the same if the design values of Table I1 are used. Compared to any standard steady-state simulator, DYNSIM needs Uotherninformation related to the process units. In steady-state simulations with a standard simulator, for a PT flash problem (for a closed system), only the temperature, the pressure, and the feed composition need to be given. For the same problem in a dynamic simulator like DYNSIM, sizing parameters and manipulative variables need to be specified instead of the temperature and pressure. Thus, in order to match the results from a steady-state PT flash simulation, the corresponding sizing parameters and manipulative variables (in this case, the external heat duty term) have to be determined. It should be noted that the sizing parameters do not affect
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 249 >
ulation can be matched. First the conditions for a reference steady state are obtained through a combination of dynamic and steady-state simulations with and without energy balance. That is, the simulation ie started without energy balance with the dynamic mode, when a condition close to the steady state is reached, the energy balance is added, and after continuing for a period of time, a switch to the steady-state mode is made. At the end of this simulation, the steady-state conditions that are obtained are stored for future use. The effect of disturbances/changes in the design variables are studied with respect to this reference steady state. That is, the simulation strategy shown in Table I for problem type ii has been followed. Parta a-c of Figure 3 show the transient responses for disturbances in manipulative variables (feed flow rate and heat duty) and the sizing parameter (valve coefficient for liquid flow). As one would expect, changes in the sizing I b variables cause an initial disturbance in the process, but after a period of time the process comes back to the same steady state. The transient response of pressure for flash unit 1shows an inverse response for a disturbance in the I I heat duty to the flash unit (see Figure 3c for +lo% change in heat duty). These results come from solutions of EVAPORATOR MIXER problems of types i and ii. Figure 3d shows the effect of changes in the volume of flash tank 2. It can be seen that Figure 2. (a) Flow sheet for the modified Cavett problem. (b) Flow the process comes back to the same reference condition. sheet for the 1-hexene process. This problem is interesting in the sense that since the flash tanks operate near the two-phase boundaries,they are very the final steady-state condition but they do affect the the sensitive to changes in the design variables even if the final transient path leading to this steady state. steady state is the same. The final steady states have been In this work, therefore, we have designed the flash units confirmed by solving problems of type iii. and the mixers such that the results from a PT flash simb LIQUID FLOW FROM FLASH NO. 4 (KMOL/HR) VAPOR FLOW FROM FLASH NO. 2 (KMOL/HR) a
PROD. 1
--Ezl-D
130.00
100.00
4
if
&
125.00
I
>
0 120.00
.I Y
>
_1
70.00
-
115.00
-
+lo$
+lox IN FEED FLOW + l o x IN HEAT DUTY
PRESSURE IN FLASH
C
2.04
G 2.00
NO.
IN HEAT DUTY
tlh IN LIQ VALVE
+ 1 h IN LIP VALVE
+ 1 h IN FEED FLOW
'
1 (ATM)
d
PRESSURE IN FLASH NO. 2 (ATM)
3
,
I
0 000
,
I
1
#
,
I
0.002
t
/
,
,,
I I I / I / ,8 1 8
0.004
8
7
I I I I I
0.006
~ " " ' I ' ' " ' " ' ' I
0.008
0.010
TIME IN HRS.
- _
t10x IN HEAT DUTY
2 + l o % IN LIC VALVE ...... + i o % I N FEED FLOW'
-
3% INCREASE IN VOLUME _ _ -10% INCREASE IN VOLUME
...... 17% INCREASE IN VOLUME,
Figure 3. (a-c) Transient responses due to changes/disturbances in design parameters and manipulative variables for the Cavett problem. (d) Transient responses due to changes in the volume of flash 2 for the modified Cavett problem.
250 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
l-Hexene Production. This process includes a plug flow reactor, an evaporator, a flash tank, and a mixer. l-Propene is converted to l-hexene in a plug flow reactor. The feed contains an equimolar mixture of l-propene and propane. The flow sheet for this process is shown in Figure 2b, and the specified values of the design variables are given in Table 11. The Soave-Redlich-Kwong (SRK) equation of state (Soave, 1972) has been used for phase equilibrium predictions. Simulations are started from a condition where l-hexene is not present in the system. The simulation strategies of Table I have been employed. First a type i problem is solved until t = t4 and then a type iii problem is solved by using known values of y, x, and p (some element of p are specified) and computing the unknown elements of d. An example of the study (problem of type iii) of the transient responses and the final steady state as a function of equipment-sizing parameters is shown in Figure 4a-c. In this figure, two different values for the area of the evaporator have been used. In this figure, the steady-state value has been obtained by solving problems of type iii and the transient responses have been obtained by solving problems of type i. Aromatic Extraction Process. Gani et al. (1990)have verified the applicability of DYNSIM with respect to simulation of an aromatic extraction process consisting of highly interlinked separation columns (liquid-liquid extraction, distillation), heat exchangers, and mixers/splitters/decanters. In this process (for flow sheet, refer to Gani et al. (1990))a mixture of paraffins and aromatics are separated by using a solvent mixture. The process starta with a liquid-liquid extraction where the solvent mixture separates the aromatics from the paraffins. The extract containing the solvent mixture, the extracted aromatics, and some paraffms is stripped of the remaining paraffins in a stripping column, and the bottom product of this column is sent to a solvent recovery column. From the solvent recovery column, the solvent is mixed with make-up solvent and recycled to the liquid-liquid extractor. The process is highly interactive because the bottom product from the solvent recovery column supplies heat to the stripping column. Also, the top product from the stripping column is split into two streams. One stream is taken out, and the other stream is mixed with the solvent going into the liquid-liquid extractor. For this paper, the paraffins are represented by heptane and the aromatics by benzene and the solvent mixture consists of sulfolane and water. Simulations have confirmed that this process is highly interlinked and small changes/disturbances in the operating condition affect the entire process. In this paper, the effects of disturbances on a selection of manipulative variables are studied and a simple solution is proposed to minimize the effect of the interactions. Some of the transient responses with and without the proposed simple solution are shown in Figure 5a-d. Run 3 corresponds to simulation without a mixer between the two distillation columns for a disturbance in the reboiler heat duty of the first distillation column. Runs 1, 2, 4, and 5 correspond to simulations with a mixer between the two distillation columns. The effect of adding the mixing tank is to reduce the interactions between the different process units. The simulation results, as shown in the Figure 5 a 4 , confirm this. These simulation results were obtained with the simulation strategies for problems of types i and ii from t = t4-t6 (see Table I). Near-Critical Extraction of 2-Propanol from Water by Isobutane. Brignole et al. (1987) have presented steady-state simulation and design for this process. In this paper, we have simulated this process with DYNSIM using
-
3
- A =
4
vi
- _
A =
...
STEADY STATE
27 M
LIQUID FLOW FROM FLASH USING VARIOUS DESIGN b
3
4000
000
0 1c
020
030
040
050
T M E (rRS:
-----A= - - A
4
MZ
= 2 7 MZ STEADY S'ATE
PRESSURE IN EVAPORATOR USING VARIOUS DESIGN C
' 00°0 6
300 200
1 1
1
-7n
000
, I , , , , , ) , , ,
I I , , , , )
0 10
__
_-
n 050
I , I I / , , , , ~ , I , , / , , , I ) , , , , , I / I
020 030 TIME (HRS)
040
A = 4 hj' A = 27 M STEADY STATE
Figure 4. (a-c) Transient responses due to changes in the area of the evaporator for the l-hexene process.
the same thermodynamic model used by Brignole et al. (1987). The process flow sheet is shown in Figure 6. We are able to obtain the same steady state as reported by Brignole et al. (1987). Since in many design and operability studies, simulations around a reference steady state are needed, we have tried to investigate further the problem-specific simulation strategies that have been proposed by Gani et al. (1990). The distillation column (2) is selected for this analysis. First we obtain a reference steady-state condition. At this condition, we fit the equilibrium-phase compositions on each plate to the
Ind. Eng. Chem. Res., Vol. 31, No. 1,1992 261 a
Run 3
,
b
,
Run 3
0 514 7 C912 N
c
5
C 910 0 908 Run 4
L
0
-
3 906
Run 5
0
0 9c4
Run 2 Run 1
0
2 0.200
0
L
c
.-C -0.190
3
- - -*-
?
v
T O 902 v
-
---_-.
Run 5
X
0 900
X 0.185
J
0 838
0.180
1
I " ' ~ 1 ' " ' I 8 ~1 ~ ~ '~" I " ' , l ~ ~ I I / ~ ~ , ' l
0 0 0 . 3 0.5 0.8 1.0 1, 3
0.0 0.3 0.5 0.8 1.0 1.3 1.5 1.8 2.0 2 . 3 2.5
,
C
1.5, 1.8 2.0 2.3 2 5
Time in mtn.
Time in min. Run 3
I
0.182 7
0.181
,, , ,
I
0.180 3
50.179 0 L
0.1 7 8
0
Run 5
E 0'177 0 0 176
+ d
Run
80.175
Run 2 Run 4
.E 0 . 1 7 4
0 172
?
0 171 0.0
0.3
0.5 0.8 1.0 1,.3 1.5 1.8 2.0 2 . 3 2.5
0,C 0.3 0.5 0.8. 1.0 1..3 1.5, 1.8 2.0 2.3 2 5
Time in min.
Time in min.
Figure 5. (a-d) Transient responses of variables of the distillation columns due to changes in the operating condition: (run 1) reboiler heat duty of column 1 reduced by 2.3%;(run 2) reboiler heat duty of column 2 reduced by 1.4%;(run 3) reboiler heat duty of column 1 reduced by 2.3% (without mixer between columns 1 and 2); (run 4) reflux rate of column 2 reduced by 1.1%;(run 5) heptane (1)-benzene (2) feed to the extractor changed from 0.25, 0.75 to 0.23, 0.77 (mole fraction). Table 111. Coefficients for Model Simplification (Equation 15)
1: 55.5 moles/h 2: 16.6 T = 332.6 K
water (1) 500 'I = 300
component water 2-propanol isobutane water 2-propanol isobutane water 2-propanol isobutane
K
Pumps are not showo. IS ' T is stream conditioner.
DIST2 1
3: 48.0
1: 1.5 mol(:s/b
section
coefficients A B 7.1399 -2920.997 9.5419 -4137.923 6.4215 -1856.208 2.7625 -1293.450 10.145 -4362.244 9.3874 -2959.660 9.0096 -3722.240 8.7793 -3810.075 2.7334 -270.7030
temp interval/ K 360-435 360-435 360-435 360-400 360-400 360-400 400-435 400-435 400-435
in Figure 7a-d. The reference steady-state profiles are shown in Figures 8a-c. The coefficients for aiand bi are given in Table III,and the column specificationsare given in Table IV. In Figure 7a, the response for the simplified model shows an inverse response. This is probably not real but due to the coefficients for eq 15. In Figure 8a, it can be seen that the differences between the actual values and the fitted values are larger at the top end of the column. The dotted line and the asterisk (Figure 8a-c) show the result from two different fitting procedures. The asterisk corresponds to one correlation for the entire column, and
252 Ind. Eng. Chem. Res., Vol. 31, No. 1,1992
-15
:I, ;>,
-20
,I
-35
0.8
1:2
1
-14:
,
, , ~, I ~, , I , , , , , , L6, 2.0 2 . 4 2.8 3.2
, ,
0.0 0.4
Full model
,
,
-16
1
C: 0
3.6
Time in min.
0
1
0.3
0.2
01
05
04
0.7
0.6
Time in min.
-0.0 c5
Full model
0
Simplified model
__c___________....._......~..
Full model
--_. Simplified model
-5.3 -5.5
~
E0
1
,
,
,
1
,
1.0
,
,
,
1
,
,
,
,
,
,
,
2.0 3.0
,
,
,
,
,
,
40
,
1
,
,
,
/
1
,
/
,
,
1
,
,
,
,
-0.15
,
5.0 6.0 7.0
8.0
~
00
,
1
,
,
~
/
10
,
1
,
~
2.0
,
,
,
,
,
3.0
,
,
,
,
~
4.0
,
/
,
1
/
5.0
,
,
,
,
1
,
,
,
,
,
,
6.0 7.0
,
,
/
,
8.0
Time in min. Time in min. Figure 7. (a,b) Comparison of transient responses (and steady states) of the vapor cornposition of isobutane from top plate &(I)) and the temperature of plate 5 (2'5) due to a disturbance of +5% in the feed flow rate. (c,d) Comparieon of the transient reaponees (and steady states) of the liquid composition of water in the bottom plate ( ~ ~ ( 1 0and ) ) the temperature of the bottom plate due to a disturbance of +5% in the reboiler heat duty. Table IV. Column Specifications for the Near-Critical Extraction Process variable value no. of plates (actual) 22 no. of plates (theoretical) 10 reflux rate/(kmol/h) 30.0 2500.0 reboiler heat duty Qb/(KJ/h) feed stream location (plate no.) state composition water 2-propanol isobutane temperature/# pressure/atm flow rate/(kmol/h)
0.023 0.250 0.727 347.0 11.0 100.0
plate geometry detaila type of plate weir length/m weir height/m active area/m2
sieve 0.53 0.038 0.366
5 liquid
the dotted line corresponds to two correlations-one correlation for above the feed plate and one for below the feed plate. These results show that even for highly nonideal mixtures it is possible to simplify the implicit algebraic equations (eq 2) and obtain satisfactory results (in terms of percentage error).
Table V. Comparison of Computing Times for the NCE Simulation Problem computational time linearized problem full model model ODE DAE step in feed 1.00 0.250 1.00 (sparse) 2.465 1.00 (dense) 1.120 step in Qb 1.00 0.278
The computing times corresponding to simulation with the fullmodel and the simplified model are given in Table V. It can be seen that nearly 75% of the computing time can be saved for this process by using eq 15. A large savings is possible because the GCEOS model is computationally more expensive than other thermodynamic models (for example, the SRK (Soave, 1972) equation of state). The computing times for the full model are compared also for the ODE mode with sparse technique and the DAE mode. The DAE mode needed more computing time than the ODE mode. This is because the procedures need very few iterations in the ODE mode (for distillation columns), while in the DAE mode, the size of the problem is nearly twice as much as that in the ODE mode.
Conclusions Design and analysis of chemical processes can be performed through process simulation. Since the behavior of the model equations varies with the problem type, the
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 253 I
440
I
420 41 0 y 4 0 0 v
390 Q
E 380 a, -I
370 360
:::
,,,,/,,,, /,,,,
\ 7 , , l
( , , /
~
, ( , / ,
I
, , ( , , ,/ , , I , , , ,
330 1
2
5
4
5
6
7
8
9
1
0
Plate No. b 1 .oo
0.90 ~ 0 . 8 0 0 > 0.70
.-+-0 0.50 0.40
2-Propanol 0.10 0.00
1
2
3
5
4
6
7
8
9
10
8
9
10
Plate No C
100
: 2-Propanol
090 1
.& 0
60
i
0 0.50
80
40
lsobutane
7
-
.
-
4
7
5
6
7
Plate No. Figure 8. (a-c) Temperature and liquid and vapor composition profilea for distillation column 2 (#) fitted with one correlation for
equation-oriented simulators like DYNSIM and that the proposed methods of solution are ideally suited for integration of simulation and design. The results therefore serve as examples of integration of process simulation (dynamic and steady state), design, and analysis. It should be added that not all of the simulation strategies mentioned in this work have been implemented in DYNSIM to work in the automatic mode. The different switching times (between dynamic and steady state, stiff and nonstiff, and ODE and DAE modes) take place at user-specified values or default values set by DYNSIM. Future work will provide an automatic switching procedure with variable times. Currently work is in progress to test the applicability of new thermodynamic models (for example, the MHV2 model of Dahl et al., 1991),which should be more accurate and faster than the GCEOS model and ideally suited for studies related to supercritical extraction (SCE) processes. The SCE process for dehydration of acetone (Cygnarowicz and Seider, 1989) is currently being studied along with several other highly complex chemical processes, for example, design/analysis of heterogeneous azeotropic distillation processes.
Nomenclature ai = coefficient for component i (eq 15) A = cross-sectional area, m2 bi = coefficient for component i (eq 15) dl = equipment sizing parameters d2 = manipulative variables E = activation energy, MJ/kmol F = flow rate, kmol/h h = enthalpy, MJ/kmol HR = heat of reaction, MJ/kmol J = Jacobian matrix KO = reaction rate constant, (kmol/m3)'-"/h for nth order reaction Ki = equilibrium constant of component i ni = molar holdup of component i, kmol p1 = explicit algebraic variables (function of y and d) p2 = explicit algebraic variables (function of x and d) p o = porosity (fraction) P = pressure, atm Pmh= minimum allowable pressure, atm (see Table 11) Q = heat duty, MJ/h Rm = valve coefficient for mixer RL = valve coefficient for liquid flow from flash tank RV = valve coefficient for vapor flow from flash tank t = time, h T = temperature, K u = heat-transfer coefficient, MJ/(m2 K-l h-' 1 ui = molar volume of component i, m3 V = volume, m3 x = implicit algebraic variables x i = mole fraction of component i y = differential/algebraic variables Subscripts i = component i L = liquid phase V = vapor phase
the entire column; (- -) fitted with two correlations for two sections of the column; (-) reference condition.
Literature Cited
inherent flexibility of the equation-oriented approach offers scope for the adoption of different simulation strategies. The simulator to be used, however, must have special features which will make it possible to solve the different simulation problems related to design and analysis. The simulation results confirm that different aspects of the design problem can be studied through
Brignole, E. A.; Andersen, P. M.; Fredemlund, Aa Supercritical fluid extraction of alcohols from water. Ind Eng. Chem. Res. 1987,26, 254-261. Cameron, I. T.; Gani, R. Adaptive Runge-Kutta algorithms for dynamic simulation. Conput. Chem. Eng. 1988,12 (71,705-717. Cavett, R. H.Application of numerical methods to the convergence of simulated processes involving recycle loops; American Petroleum Institute Preprint No. 04-63;American Petroleum Institute: Washington, DC, 1963.
254
Ind. Eng. Chem. Res. 1992,31, 254-262
Cygnarowin, M. L.; Seider, W. D. Effect of retrograde solubility on the design optimization of supercritical extraction processes. Znd. Eng. Chem. Res. 1989,2%, 1497-1503. Dahl, S.; Rasmussen, P.; Fredenslund, Aa. The MHV2 model: A UNIFAC-based equation of state model for prediction of gas solubility and vapor-liquid equilibria at low and high pressures. Znd. Eng. Chem. Res. 1991,30,1936-1944. Duff, I. S. A set of FORTRAN subroutines for sparse unsymmetric linear equations. U.K.Atomic Energy Authority, Harewell Laboratory, [Report] AERE-R; AERE-R8730; Her Magisty’s Stationary Office, HMSO London, 1979. Fredenslund, Aa.; Gmehling, J.; Rasmussen, P. Vapor-Liquid Equilibria Using UNIFAC. Elsivier: Amsterdam, 1977. Gani, R.; Cameron, I. T. Extension of dynamic models of distillation columns to steady state sirnulation. Comput. Chem. Eng. 1989, 13 (3), 271-280.
Gani, R.; Perregaard, J.; Johanaen, H. Simulation strategies for design and analysis of complex chemical processes. Trans. Znst. Chem. Eng. 1990,68 (PartA), 407-417. Hillestad, M.; Hertzberg, T. Dynamic simulation of chemical engineering systems by the sequential modular approach. Comput. Chern. Eng. 1986,10,377-388. Pantelides, G. C.; Gritsis, D.; Morison, K. R.; Sargent, R. W. H. The mathematical modelling of transient systems with differentialalgebraic equations. Comput. Chem. Eng. 1988,12 (5), 449-454.
Perregaard, J.; Pedersen, B. S.; Gani, R. Steady state and dynamic simulation of complex chemical processes. Presented at the 4th International Symposium on Process Systems Engineering, Montebello, Canada, Aug 5-9, 1991. SkjoldJlargensen, S. Group contribution equation of state (GCEOS): A predictive model for phase equilibrium computations over wide ranges of temperature and pressures up to 30 MPa. Znd. Eng. Chem. Res. 1988,27,110-118. Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 1972,27,1197-1203. Smensen, E. L. Dynamic simulation of chemical processes with accurate thermodynamic models. Ph.D. Thesis, The Technical University of Denmark, September, 1991. Ssrensen, E. L.; Johansen, H.; Gani, R.; Fredenslund, Aa. A dynamic simulator for design and analysis of chemical processes. In Process Technology Proceedings; Bussemaker, H. Th.,Iedema, P. D., E&.; Elsivier: Amsterdam, 1990; Vol. 9, pp 13-18. Tiegs, D.; Rmmwen, P.; Fredenslund, Aa. Vapor-liquid equilibria by UNIFAC group contribution. 4. Revision and extension. Znd. Eng. Chem. Res. 1987,26, 159-161. Received for review March 11, 1991 Revised manuscript received September 3, 1991 Accepted September 12,1991
Gross Error Detection in Serially Correlated Process Data. 2. Dynamic Systems Chen-Shan Kao, Ajit C. Tamhane, and Richard S. H. Mah* Robert R. McCormick School of Engineering and Applied Science, Northwestern University, Evanst on, I1linois 60208- 3120
The importance of gross error detection and identification is widely acknowledged, but so far relatively few publications in the chemical engineering literature have addressed this problem for stochastic dynamic systems. Also, statistical process control (SPC) chart techniques have been developed and used extensively to detect shifts (a change in steady state) in the process means but have been hitherto ignored in applications involving detection of gross errors in chemical process data. In this paper a composite test procedure for detecting and identifying gross errors is proposed: we first make use of the CUSUM test (which underlies the CUSUM chart used in quality control applications) to detect the presence of gross errors. We then apply the generalized likelihood ratio (GLR) method to identify as well as to estimate the magnitudes of the gross errors. We concentrate our treatment on discrete-time, linear dynamic systems operated around a nominal steady state and corrupted by process and measurement noises, which are assumed to be Gaussian. These noises can be white or serially correlated. Computer simulation studies and some analytical results are presented. Methods for gross error detection and identification in steady-state chemical processes have been studied extensively in the past two decades (Tamhane and Mah, 1985; Mah, 1990). But in reality, even for steady-state operations, process conditions vary continually about a nominal steady state. A truly steady state is almost never attained. A dynamic system will be a better representation of a real process. A mathematical model which describes the dynamic behavior of the process under normal operating conditions based on physical information and statistical data is a very powerful tool in gross error detection and in process performance evaluation. Detection of gross errors in a dynamic system is more difficult than in a steady-state system since the effects of gross errors at any time depend on the dynamic behavior of the process, and thus may change with time. In a previous paper (Kao et al., 1990),we investigated the effect of serial correlations of process measurements
* To whom correspondence should be addressed. 0888-5885 f 92/2631-O254$03.Oo f 0
in the steady-state case. We showed the extreme sensitivity of the measurement test (MT) for gross error detection to serial correlations. We also proposed two methods which work for serially correlated data. The first one involves suitably adjusting the variance of the test statistics. The second, prewhitening, involves filtering out the serial correlations and then applying the desired test based on the independence assumption. In dynamic systems, serial correlations can arise from process and measurement noises. Therefore it is of interest to extend the previously developed methods of gross error detection (particularly the prewhitening method) to the dynamic case. This paper represents a major step in this direction. Only a few publications in the chemical engineering literature have addressed gross error detection for dynamic systems (Behgham and Lees, 1977;Watanabe and Himmelblau, 1982;Narasimhan and Mah,1988;Gertler and Luo, 1989). We begin by considering a discrete-time, linear dynamic system operated around a nominal steady state and cor@ 1992 American Chemical Society
