Methodology for the Design and Analysis of ... - ACS Publications

perform the design, control, and analysis of chemical processes, ..... Figure 5. Molar selectivity (SG/H) and recycle flow rate (f7), each as a ...
0 downloads 0 Views 422KB Size
8084

Ind. Eng. Chem. Res. 2007, 46, 8084-8100

Methodology for the Design and Analysis of Reaction-Separation Systems with Recycle. 2. Design and Control Integration Edgar Ramı´rez*,† and Rafiqul Gani CAPEC, Department of Chemical Engineering, Technical UniVersity of Denmark, DK-2800 Lyngby, Denmark

The simultaneous solution of the design and control problems of chemical processes represents a complex task, given that several aspects should be fulfilled (e.g., although a design can be considered to be feasible, when attempting to control it might not be possible, or the control configuration could be quite complex). The interaction between design and process variables also must be considered while attempting to outline a control configuration. The second part of this communication shows how, through a model-based methodology, it is possible to solve the design and control problems in a systematic way, leading to a full understanding of a chemical process. The Tennessee-Eastman (TE) problem has been selected to highlight the application of the model-based methodology in a plantwide problem. 1. Introduction As discussed in the first part of this series,1 a hierarchical three-stage model-based Methodology has been developed to perform the design, control, and analysis of chemical processes, showing a reaction-separation with recycle (RSR) scheme. The Methodology makes use of the potential that a model-based analysis of the process offers, with the objective of revelaing the main interactions that exist between the different design and process variables, so that the design and control integration can be performed systematically. The design perspective of the model-based Methodology was presented through a case study (the production of ethylbenzene), i.e., the analysis of regions where the process showed a desirable (static) behavior were considered only, as well as different scenarios where, for example, a “snowball effect” could be present and the values for the corresponding design/process variables are involved. One of our main concerns in this paper is how to determine an appropriate set of design (manipulated) and process (controlled) variables that lead to a suitable control structure sketch. For instance, Larsson and Skogestad2 defined the control structure design as the structural decisions involved in control systems design, including (see also Nishida3) the selection of (i) controlled variables cs (“outputs”), (ii) manipulated variables m (“inputs”), (iii) measurements ν, (iv) the control configuration (a structure interconnecting measurements/setpoints and manipulated variables (i.e., the structure of the controller K which interconnects the variables cs and ν (controller inputs) with the variables m)), and (v) controller type (control law specification, e.g., proportional-integral differential (PID), decoupler, linear quadratic Gaussian (LQG), etc.). On the other hand, defining the control structure design problem mathematically is difficult,2 and an alternative method involves the use of heuristic-based approaches, relying on experience and process understanding. In this respect, Fisher et al. presented a hierarchical heuristicbased approach that addresses the controllability,4 operability,5 and selection of a set of controlled variables6 for a complete * To whom correspondence should be addressed. Tel.: (+52) 5557296000, ext. 55394. Fax: (+52) 5555862728. E-mail: eramirezj@ ipn.mx. † Present address: Department of Chemical Petroleum Engineering, ESIQIE, IPN, C.P. 07738, Mexico City, Mexico.

chemical plant; several guidelines are given for the different tasks to be performed. The book by Luyben et al.7 recommends a procedure for the design of a control structure; however, they do not establish guidelines to systematically determine which variables must be manipulated and controlled. Recently, with the objective of circumventing this key issue, Skogestad8 presented a systematic procedure for the plantwide control structure design where the selection of the controlled and manipulated variables is based on process insights as well as heuristics/process knowledge. The use of “generic” models is also emphasized, where the structural part is correct, even though all the parameters may not match the true plant in question. A first-principles model, based on material and energy balances, is recommended, given that a good control structure is insensitive to parameter changes. Kiss et al.9 based their approach on a nonlinear analysis of the plant to identify regions in the parameter space where the system has a desired behavior. Thus, given the interaction between the design and control of recycle systems, the type of behavior exhibited by the entire plant is influenced by the plantwide control structure initially selected. On the other hand, optimization-based methods have been applied to address the plantwide control problem where the main objective has been to synthesize a design and a control structure that satisfies a set of economic constraints and operability issues. Luyben and Floudas10,11 proposed a method to consider steadystate economics and open-loop controllability objectives within the mathematical programming framework of process synthesis simultaneously. In the Mohideen et al.12 framework, the process system is modeled using dynamic mathematical models where variations include both uncertain parameters and time-varying disturbances while control structure selection and controller design is considered to be part of the resulting stochastic mixedinteger optimal control mathematical formulation (MIDO). The disadvantage of this approach is that the resulting mixed-integer nonlinear programs (MINLPs) are very large, even for relatively small-scale differential-algebraic equation (DAE) systems, which potentially prohibits the solution of large, realistically modeled systems. Mo¨nnigmann and Marquardt13 presented an approach to incorporate constraints in the system dynamics into the optimization-based design of nonlinear systems. In this approach, boundaries in the parameter space of the nonlinear system, where desired dynamics characteristics of the system are lost, are considered. Parametric robustness is guaranteed by staying off the boundaries at a user-specified distance in the

10.1021/ie061328p CCC: $37.00 © 2007 American Chemical Society Published on Web 10/20/2007

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8085

Figure 1. Tennessee-Eastman (TE) problem flowsheet.

Figure 2. Simplified flowsheet for the TE problem.

parameter space, which ensures that the desired dynamic characteristics will be met, despite parametric uncertainty or parameter drift. The aforementioned discussion gives us the opportunity to present the complete Methodology, with respect to the design and control integration. It is also important to remark that the present work is not a control study as such; however, through the application of the Methodology, sufficient design/control information is provided to support further control-related studies. Even though we argue that we will look for an “optimal solution” (within the context of the Methodology), the resulting design and control integration might be a local solution to the problem at hand, based on the design target objectives. On the other hand, an interesting way to make the proposal of a control structure is performed, based on the corresponding impact that design and process variables have on the process, which helps to ensure the design targets of the process. In this view, this work is divided in the following parts. First, a brief summary of the model-based Methodology (presented earlier1) will be given, which will lead to its application in a particular case study (the TE problem). Finally, we will draw some conclusions. 2. The Model-Based Methodology The model-based Methodology is comprised of three successive stages, where special-purpose models are used for

parametric sensitivity analyses, especially in the first two stages, and validation in the third stage. Because the integration of the design of the process and its control structure involves solving these two problems in the early stages of the design process simultaneously, it is necessary to use appropriate models that represent specific features (operation, behavior) of the process and the variables that may affect them. In stage 1, some of the design variables are combined together into coupled parameters or dimensionless numbers (such as the Damko¨hler number, Da). Therefore, the models in stage 1 are based on these coupled variables, and the objective is to determine the values of these variables where the process has the optimal behavior. Therefore, for different values of the design variables (e.g., flow rates, temperatures, etc.), the model equations are solved to obtain values of the process variables (compositions, pressures, yields, etc.), which are then used to estimate the variables that determine the process behavior (performance). This helps to decide target values for the coupled variables for the design variables (corresponding to the desired process performance), the operational limits (in terms of the process variables), and desired (target) of the process behavior. All these variables affect the design of the process, as well as the design of the control structure. In stage 2, the coupled parameters or dimensionless numbers are disassembled, thereby generating a new model, which obviously is more complex than the model in stage 1. The objective here is to determine the design variables that can match the process behavior (defined by the process variables) and the target values of the stage 1 design variables. If the stage 1 design targets cannot be matched, new targets are set until a match is found. The results from the model analysis determine the operational feasibility of the process, as well as the search space for the design-process variables. In stage 2, the overall design targets are set and the values for a minimum set of design and process variables that matches this target is found. This means that all of the decisions that are related to the design-control structure have been made in the early stages of the design process. The process can be described in terms of the performance it will give, the values of the design variables that can achieve it, and the values of the process variables that will be attained by the process. From a control structure point of view, these process variables will need to be controlled by manipulating the identified design variables. Also, these calculated values correspond to the setpoints for the manipulated

8086

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Figure 3. Conversion of A and selectivity, each as a function of the Damko¨hler number (Da): (a) xA vs Da and (b) SG/H vs Da.

Figure 4. Conversion of A and molar selectivity, each as a function of the Da number: (a) xA vs Da and (b) SG/H vs Da.

and controlled variables. The models in stages 1 and 2 are simple, because they only include the features of the process that must be studied. In stage 3, a rigorous model of the process is used to validate the design decisions and, if necessary, improve the design. The results from stages 1 and 2 are used as initial estimates; therefore, the robustness and efficiency of the rigorous model-based solver is improved. Typically, steady state as well as dynamic simulation models can be used at this stage, whereas in stages 1 and 2, mainly steadystate models are used. On the other hand, process disturbances can be chosen based on actual process information, as well as from the resulting (static and dynamic) analysis, which returns their range of variation, and also from the analysis of the degrees of freedom (DOFs). Hence, given all this information, the corresponding control proposal will be derived. 3. The Tennessee-Eastman Problem The Tennessee-Eastman (TE) problem first appeared at an American Institute of Chemical Engineers (AIChE) meeting in 199014 and at the Chemical Process Control Conference in 1991; since then, it has led to many publications and much research that involves it. This is a good example of what plantwide control involves, including nonlinearities, not only in the reaction section but also due to the interconnections. The perspective to be taken in this study is related to the analysis of its behavior from a design viewpoint and, from there, attempt to sketch a proposal for the control structure. It is important to remark that the objective of this work is not to perform a control study of the process as such, but to provide sufficient information to outline an appropriate control structure for the problem. 3.1. Process Description. The process produces two products from four reactants. Also present are an inert and a byproduct,

making a total of eight components: A, B, C, D, E, F, G, and H. The original reaction scheme is given as follows:

A(g) + C(g) + D(g) f G(l) (Product 1) A(g) + C(g) + E(g) f H(l) (Product 2) A(g) + E(g) f F(l) (byproduct) 3D(g) f 2F(l) (byproduct) All the reactions are irreversible and exothermic. The reaction rates are a function of temperature and represented through an Arrhenius expression. The reaction to produce G has a higher activation energy, resulting in higher sensitivity to temperature. Also, the reactions are approximately first-order, with respect to (wrt) the reactant concentrations. The process has five major unit operations: the reactor, the product condenser, a vapor-liquid separator, a recycle compressor, and a product stripper. Figure 1 shows the flow diagram of the process. The gaseous reactants are fed into the reactor, where they react to form liquid products. The gas-phase reactions are catalyzed by a nonvolatile catalyst dissolved in the liquid phase. A continuously stirred tank reactor (CSTR) has an internal cooling bundle, which is used to remove the heat of reaction. The products leave the reactor as vapors, along with the unreacted reactants; the catalyst remains in the reactor. The reactor product stream passes through a cooler, to condense the products, and, from there, flow into a vaporliquid separator. Noncondensed components are recycled back through a centrifugal compressor to the reactor feed. Condensed components move to a product stripping column to remove remaining reactants by stripping with feed stream 4. Products G and H exit the stripper base and are separated in a downstream refining section which is not included in the problem. The inert

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8087

and byproduct are primarily purged from the system as a vapor from the vapor-liquid separator.

0)1-

The reaction kinetics are as follows:

(

)

σ - RA,7 yA,3f3 - Da(Rˆ 1 + Rˆ 2 + νRˆ 3) (4a) σ

)

0 ) yC,1f1 -

) )

R2 ) 2Vv,r exp 10.27 -

19500 p 1.15pC,r0.370pE,r1.00 (2) RgTr A,r

0 ) yD,1f1 -

R3 ) 3Vv,r exp 59.50 -

59500 p (0.77pD,r + pE,r) (3) RgTr A,r

0 ) yE,1f1 -

where the byproduct reactions have been added to produce eq 3. Therefore, the modified scheme of reaction is

A(g) + C(g) + D(g) f G(l) (Product 1) NC

A(g) + C(g) + E(g) f H(l) (Product 2) νA(g) + D(g) + νE(g) f F(l) (Byproduct)

0)

Step 1.1: Flowsheet Simplification. The first step requires simplification of the process flowsheet, in terms of the main operations (which are mixing, reaction, and separation), which provides the basis for the mass and energy balance models (see Figure 2; nomenclature for this stage will also be referenced to this figure). Step 1.2: Reaction Analysis. The following main assumptions are made:

(4b)

σ - RC,7 yC,3f3 - Da(Rˆ 1 + Rˆ 2) σ

(4c)

σ - RD,7 yD,3f3 - Da(Rˆ 1 + Rˆ 3) σ

(4d)

σ - RE,7 yE,3f3 - Da(Rˆ 2 + νRˆ 3) σ

(4e)

(A2) The feed flow rate of component A (FA,1) is taken as a variable of the reference for the dimensionless model. (A3) No recycle of products or byproducts back to the reactor is considered. (A4) The important products are G and H, so their yields are followed, as well as the selectivity of G, with respect to H (SG/H). Step 1.3: Thermal Effects. Because the thermal effects are particularly important, the model-based Methodology recommends that we move to step 1.5. Step 1.4: Model Development Based on Mass Balances Only. Not necessary. Step 1.5: Isothermal Assumption. Although the thermal effects are important, it has been considered useful to first perform an analysis with only the mass balance model. This is done with the objective of determining the effect of, mainly, the recycle stream. (This assumption would imply perfect control of temperature and pressure in the reactor.) Step 1.6: Model Development. Equations 4 describe the component mass balances at steady state for the simplified system; thus,

( ( (

) ) )

0 ) -yF,3f3 + Da(Rˆ 3)

(4f)

0 ) -yG,3f3 + Da(Rˆ 1)

(4g)

0 ) -yH,3f3 + Da(Rˆ 2)

(4h)

( ) Ri,7 σ

- 1 yi,3f3 + Da(Rˆ 1 + Rˆ 2 + Rˆ 3)

(4i)

Rˆ 1 ) yA,31.08 yC,30.311 yD,30.874

(5a)

Rˆ 2 )

2 exp[c1,2 - c2,2/(RgTr)] 2.52 yA,31.15 yC,30.370 yE,31.00 P 1 exp[c1,1 - c2,1/(RgTr)] r

(5b)

Rˆ 3 )

3 exp[c1,3 - c2,3/(RgTr)] 2 P y (0.77yD,3 + yE,3) 1 exp[c1,1 - c2,1/(RgTr)] r A,3

(5c)

and Da is the plant Damko¨hler number, which represents the relationship between the rate of reaction, wrt the rate of fresh feed of A:

(A0)Constant pressure within the reactor, thus, pi,r ) yi,rPr. (A1) Reaction 1 (eq 1) is taken as reference, because it has the highest sensitivity to temperature.15

)

where

where ν ) 1/3. 3.2. Stage 1. Model-Based Analysis by Simple Models. This analysis involves several steps.

NC

yi,1f1 + ∑ ∑ i)1 i)1

(

σ - RB,7 yB,3f3 σ

0 ) yB,1f1 -

42600 R1 ) 1Vv,r exp 44.06 p 1.08pC,r0.311pD,r0.874 (1) RgTr A,r

( (

(

Da )

1 exp[c1,1 - c2,1/(RgTr)]Pr2.265Vr FA,1

(6)

Step 1.6.1: Solution and Analysis of the Model Equations. Equations 7 (or, alternatively, eqs 4) represent a system of eight nonlinear algebraic equations (NLAEs), wrt the mole fractions in the reactor (yi,3) plus the overall mass balance (f3), together with the dimensionless rate of reaction Rˆ k (eqs 5); hence, Ne ) 12. To classify the variables used in the aforementioned mass balance model, a vectorial compact form of eqs 7 is given below:

0 ) f[x,u,p]

(7)

where

x ) [f3, yi,3, Rˆ k]T

(for i ) A, ..., H; k ) 1, ..., NR)

u ) [Da, σ, Ri,7]T

(for i ) A, ..., E)

(8) (9)

p ) [F1, Tr, Pr, Rg, yi,1, k, c1,k, c2,k]T (for i ) B, ... E; k ) 1, ..., NR) (10) where vector x contains the representative process variables to be solved; vector u contains the main design variables (namely, the Damko¨hler number (Da) and the recovery (Ri,7) and purge factors (σ)); and vector p represents the main reaction parameters, pressure, and temperature in the reactor, mole

8088

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

reactants enhances the selectivity of the product(s) of interest, but can also lead to operational restrictions (e.g., the snowball effect), and (ii) low Da values (Da < 50), which can be interpreted as high FA,1 (see eq 6), result in high selectivity values, because reactant A does not act as a limiting reactant, so that the overall effect is an increase in the rates of reaction. Step 1.7: Addition of Energy Balances. The corresponding energy balances for each operation are given as follows. Mixing Zone. For the simplified flowsheet, the mixing zone considers the feed stream (F1) and the recycle stream (F7) as inlets and (F2) as an outlet, thus giving the following equation for energy balance: H

0) Figure 5. Molar selectivity (SG/H) and recycle flow rate (f7), each as a function of the receovery of componetn C (RC,7).

fractions, and flow rate at the inlet stream, as well as preexponential factors and activation energies. Based on this, Nu ) 36. Therefore, the DOF is given as NDOF ) Nu - Ne ) 24. Hence, by specifying the parameters vector p and the design variables vector u for a specific set of conditions of operation, the system can be solved. The model-based process analysis is performed considering a constant pressure within the reactor (Pr ) 2.806 MPa and Tr ) 393.55 K) (the base case scenario, see Table A2 in the Appendix for process specifications). The simulation results obtained by solving eqs 4 for different sets of values of vector u are analyzed in Figures 3-5. In Figure 3a (where Ri,7 ) 0, i * A; σ ) 1), the sensitivity of the conversion of A (xA), with respect to Da, for different recovery values (RA,7) is shown. As can be seen, xA decreases because of a larger amount of A that must be reacted (due to recycle) within the same reactor volume (assuming Vr is kept constant). On the other hand, it is important to highlight the combined effect of Da and the recovery factor RA,7. That is, in Figure 3a, a steeper slope is observed in the conversion profile as RA,7 decreases, especially at lower Da values (Da < 50). This effect indicates that a faster rate of reaction is most likely occurring, because of the stoichiometric consumption of reactant A. Conversely, larger RA,7 values mean that a higher residence time in the reactor is required for A. Now, in Figure 3b, the molar selectivity in the reactor, which is defined as SG/H ) yG,3/yH,3, also is favored at low Da (Da < 10) and high RA,7 values. A value of Da > 200 has no effect on the selectivity, because it can be assumed that equilibrium has been reached. Thus, it can be argued that (i) high selectivity values are obtained at low Da values, because product G is produced faster than product H, and (ii) high recovery values increase the selectivity, because of a higher concentration of A, which makes the rates of reaction become slower (see eqs 5a and 5b), thereby causing the ratio G/H to increase. Figures 4a and 4b (where Ri,7 ) 0, i * C; σ ) 1) confirm the trend observed, with respect to RA,7. However, although it would be preferred to operate at high SG/H values, Figure 5 shows that the recycle flow rate (f7) increases rapidly as RC,7 f 1, thereby increasing the likelihood of a “snowball effect”. This is, in fact, a characteristic of the TE problem, which operates at high recovery values for unreacted components. From the aforementioned brief analysis, which assumed isothermal conditions, it was observed that (i) the recycle of

yi,jCpvi )(Tj - T2) ∑ ∑ j)1,7 i)A Fj(

(11)

Reactor. The reactor operation assumes a non-adiabatic nonisothermal condition. Therefore, heat exchange with a cooling media (water) is included in the energy balance by considering a heat balance on the reactor side: H

0 ) F3(



NR

yi,3Cpvi )(T2 - T3) -

i)A

∑ ∆Hr,kRk - Q˙ r

(12)

k)1

where H

∑ Vi,kCpvi (T3 - Tref) i)A

∆Hr,k ) HoFk +

(13)

with the heat balance on the cooling media side being given by

Q˙ r ) mcwCpcw(Tcw,r,out - Tcw,r,in)

(

Q˙ r ) UAr

∆T1 - ∆T2

ln(∆T1/∆T2)

(14)

)

(15)

∆T1 ) T3 - Tcw,r,in

(16a)

∆T2 ) T3 - Tcw,r,out

(16b)

Separation Section and Purge Zone. In this step, the reaction effects are mainly analyzed. Therefore, it is assumed that there will be no change in the outlet stream temperatures from the splitter as well as the purge; that is,

T4 ) T5 ) T3

(17)

T7 ) T6 ) T5

(18)

The dimensionless energy balance equations take, as variables of reference, the feed (inlet) to the process of component A (FA,1), because component A is considered to be the limiting reactant, and also its heat capacity in the vapor phase (CpvA). Similarly, the temperature of the inlet to the reactor (T2) is taken as a reference value, as in the definition of θj ) (Tj - T2)/T2, which is used later in eqs 19 and 20. Therefore, the dimensionless energy balance equations take the following form: Mixing Zone Model. Rearranging eq 11 for T7 ) T3 and using the definition of θj gives E

0)

E

∑ (fi,1C˜ pi)θ1 + i)A ∑ i)A

(

)

Ri,7 fi,3C ˜ pi θ3 σ

(19)

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8089

Figure 6. Conversion of A and selectivity, each as a function of Da: (a) xA vs Da and (b) SG/H vs Da.

Reactor Model. Rearranging eq 12 in terms of dimensionless variables gives E

0 ) -f3(



[∑ ( ) ] [ ] NR

C ˜ pi)θ3 + Da

i)A

k)1

C ˜ pRk exp

γk θ3

K* ˆ k (θ3 kR 1 + θ3 ∆θ1 - ∆θ2 θref) - δc (20) ln(∆θ1/∆θ2)

where ∆θ1 ) θ3 - θcw,r,in and ∆θ2 ) θ3 - θcw,r,out. The following dimensionless variables have been used in eqs 19 and 20:

1 exp(c1,1 - γ1)Pr2.265Vr Da ) FA,1 K*2 )

K*3 )

Figure 7. Dimensionless θ3 and θcw,r,out values, each as a function of Da and RA,7.

where

2 exp(c1,2 - γ2)Pr2.52 1 exp(c1,1 - γ1)Pr2.265 3 exp(c1,3 - γ3)Pr2

c2,1 RgT2

γ2 )

c2,2 RgT2

γ3 )

c2,3 RgT2

C ˜ pi ) δc )

Cpvi CpvA

UAr FA,1CpvA

mcpcw )

mcw,rCpcw FA,1CpvA

Step 1.7.1: Solution of the Equations. Equations 19 and 20, together with eqs 4, form a NLAE system that must be solved for the new set of process variables x. In the same way as that described in step 1.6, the combined mass and energy balance model is now represented in compact form as

0 ) f[x,u,p]

(22)

u ) [uMBal,δc,mcpcw]T

(23)

p ) [pMBal,T1,T2,Tcw,in,Pr,HoFk,Cpvi ]T (for i ) A, ..., H; k ) 1, 2, 3) (24)

1 exp(c1,1 - γ1)Pr2.265 γ1 )

x ) [xMBal,θ3,θcw,out]T

(21)

Here, the subscript “MBal” refers to the variables in step 1.6.1. Therefore, vector x now also includes the temperature in the reactor (θ3) and the outlet temperature of the jacket (θcw,out). Similarly, vector u adds, as design variables, a heat elimination capacity variable (δc) and cooling media variable (mcpcw); vector p represents the reaction parameters, heat capacities, etc. The total model has now 59 variables, with 19 equations (between state and algebraic equations); thus, the DOF is given as NDOF ) 40. Therefore, by specifying the design variable vector u plus the parameter vector p, the model can be solved. Table A4 in the Appendix gives the values for some of the dimensionless variables that were kept constant in the following analysis. Step 1.7.2: Model Analysis. The simulation results obtained by solving the mass balance model for different sets of values for the design variables (u) are analyzed using Figures 6-10. Da and RecoVery Factor Analyses. The results from the mass and energy balance model highlighted in Figure 6 (Ri,7 ) 0.98, i * A; σ ) 1) confirms the trends, with respect to RA,7, as observed previously with the mass balance model (see Figure 3) and, thus, the insights from step 1.6. On the other hand, because the energy balance is included, it is expected that the rates of reaction increase, which is also understood with an increase in the dimensionless variables θ3 and θcw,r,out (as a function of Da; see Figure 7). Figure 8 shows the sensitivity of the process, with respect to the recovery of C, RC,7 (at Da ) 1.0497, corresponding to the

8090

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Figure 8. Selectivity and yields of G and H, each as a function of RC,7 (Da ) 1.049): (a) molar selectivity SG/H and (b) reactor heat duty and recycle flow rate.

Figure 9. Conversion of A and yields of G and H, each as a function of fC .

base-case scenario). Note that component C is one of the main reactants and it has the largest fresh feed. Therefore, in Figure 8a, the recycle flow rate f7 increases considerably as RC,7 approaches complete recovery. Now, if the yield of component i is defined as the amount of i produced, with respect to the amount of the limiting reactant A reacted, we get

Yi )

fi,3 fA,2 - fA,3

(25)

The yields of components G and H (calculated through eq 25) are shown in Figure 8b. Here, the production of G (YG) is favored as the recovery of C increases, although component C is present in the two main reactions (note that the second reaction involves the production of H). The aforementioned analysis clearly shows that high sensitivity is expected in the main performance variables (SG/H, xA), because of the policy of high reactant recovery, thereby leading to the possibility of a snowball effect. On the other hand, the Da analysis has shown that higher selectivity and conversions are found at low Da values, as a result of the higher rates of reaction. Analysis of the Inlet Flow-Rate Ratio. Changes in the inlet flow rate might be due to an increase in the plant’s throughput or may be due to disturbances. Therefore, Figure 9 presents the simulated responses in the conversion of A and the yields of G and H, as a function of the feed ratio of component C (represented by fC(FC,1/FA,1)). The base-case scenario involves a value of fC ) 0.99 (Ri,7 ) 0.98). Note that the highest conversions of A are achieved at fC ≈ 0.99. By changing the feed ratio, it is possible to increase the selectivity, because more reactants are fed. However, an increase in Q˙ r should be expected,

given the exothermal nature of the reaction. This is confirmed through the plotted results shown in Figure 10b. On the other hand, an interesting effect is observed: if the ratio fC is reduced (from the base-case scenario), the selectivity reaches a minimum at fC ≈ 0.9, and then subsequently increases (although leading to a lower throughput). In this way, for a given selectivity SG/H, two different fC ratios can be found, as well as two different production rates. Values of fC that are outside of the range presented would be of no interest, even if we could achieve higher selectivities; Figure 10b shows that Qr and f7 will also increase, representing operational/control problems, because of the amount of energy to be removed from the reactor and the ability of the plant to handle high recycle flow rates. Step 1.9: Storage of Information. In this stage, the effect of the size of the recycle stream on the performance of the RSR system (e.g., Figure 8) has been verified. Also, how this effect is related to the Da value needed to achieve the given performance criteria (e.g., Figures 3, 8a, and 8b) has been verified. The feed ratio has been determined to be an important (design) variable. The feed ratio can be used as an indicator of the behavior of the process under disturbance or, under a shift in the operating policy (see operational window in Figure 9 and 10). 3.3. Stage 2. Model-Based Analysis by Detailed Model Analysis. This modeling stage involves the following steps. Step 2.1: Multiple Reaction Analysis. There are three parallel reactions in this process. The candidate AR area (see Figure 11) is presented only for the CSTR, because the operation involves only this equipment. The candidate AR area has been constructed under the assumption of perfect temperature control (isothermal conditions). The CSTR trajectory, with respect to SG/H, is not convex, which is a necessary condition for the AR area,16 thus requiring to add a bypass to the reactor. However, this modification has not been attempted, because the intention of the analysis is not to change the design of the TE processs. However, the results do provide directions to obtain an improved design that satisfies the criteria of the AR area. Step 2.2: Driving Force Analysis. Because of the proprietary nature of the process, sufficient information available about the thermodynamic properties of the components, or the identities of them, is not available. On the other hand, no distillationbased separation is used in this system; therefore, driVing force diagrams have not been developed. Step 2.3: Assumptions Removal. All the assumptions considered in stage 1 (A0 - A3) are removed, e.g., the pressure within the reactor is no longer a constant, and the recycle of products and byproducts will be included in the detailed model. Step 2.4: Consideration of All Phenomena. Here, all components are considered, including the mixing zone, reactor, product separator, compressor and purge, and stripper.

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8091

Figure 10. Selectivity and reactor heat duty, each as a function of fC (Da ) 1.049): (a) molar selectivity SG/H and (b) reactor heat duty and recycle flow rate.

Figure 11. Candidate attainable region (AR) for a continuously stirred tanl reactor (CSTR) for the TE-Problem.

Table 1. Manipulated Variables for Operation Modes input

base case

mode 1 (50/50)

mode 2 (10/90)

F1 F2 F3 F4 F8 F9 F10 F11

11.200 kmol/h 114.500 kmol/h 98.000 kmol/h 417.500 kmol/h 1201.500 kmol/h 15.100 kmol/h 259.500 kmol/h 211.885 kmol/h

11.991 kmol/h 114.314 kmol/h 96.471 kmol/h 413.782 kmol/h 1441.021 kmol/h 9.497 kmol/h 253.563 kmol/h 210.885 kmol/h

13.848 kmol/h 22.948 kmol/h 174.679 kmol/h 383.109 kmol/h 1419.501 kmol/h 16.164 kmol/h 243.825 kmol/h 194.638 kmol/h

mcw,r mcw,s

93.70 m3/h 49.37 m3/h

75.40 m3/h 52.02 m3/h

62.95 m3/h 47.40 m3/h

mcw,str

230.31 kg/h

4.74 kg/h

4.90 kg/h

Mixing Zone. All components are in the vapor phase; therefore, the model includes the mass and energy balance equations, plus two algebraic equations for the pressure and composition of the mixing zone. Reactor. The reactor contains liquid and vapor phases, which are assumed to be in equilibrium. There is a significant holdup of liquid of G and H in the reactor, but there is no liquid effluent. To control liquid accumulation in the reactor, one must balance production (by reaction) against vaporization and removal in the gaseous effluent (stream 7).15 The feed and the product streams of the reactor are present in the vapor phase. Product Separator. The deviation from ideality in the vaporliquid equilibrium is described by a constant activity coefficient for each (condensable) component. Generally, where liquid and vapor coexist, the accumulation of the less-volatile components (D through H) in the vapor is neglected. In other words, these components have an equilibrium partial pressure, but the corresponding number of moles in the vapor is neglected when

computing the molar holdups. Otherwise, an iterative flash calculation would be needed.17 Compressor and Purge. Ricker and Lee17 introduced the recycle stream as an independent variable. This avoids modeling the compressor in detail and the stream flow rate is no longer dependent on the valve position in the compressor recycle. Similar to the compressor simplification, the purge stream (F9) is an algebraic variable. Stripper. There is little information available about the stripper. The unit is characterized by a split fraction model, such as that described by Ricker and Lee.17 Given that the energy balance is added, split fractions (Φi) are modeled as third-degree polynomials in temperature. Step 2.5: All Process Equipment. The separation section includes, as mentioned previously, a stripper and a product separator, as well as a compressor. Step 2.6: Model Development. The model equations are based on the model by Jockenho¨vel et al.;18 however, in this part of the analysis, they are presented on a steady-state basis. Mixing zone. Within the mixing zone, all feed streams and the recycle stream are mixed and fed into the reactor. Molar balances for components A-H:

0 ) yi,1F1 + yi,2F2 + yi,3F3 + yi,5F5 + yi,8F8 - yi,6F6 (for i ) A, ..., H) (26) Energy balance for the mixing zone: H

0)

∑ Fj(i)A ∑ yi,jCpvi )(Tj - Tm) j)1, 2, 3, 5, 8

(27)

Pressure and concentrations:

Ni,m

yi,6 )

(28)

H

∑ Nj,m

j)A

R gT m Ni,m Vm i)A H

Pm )



(29)

Reactor. In the reactor model, the influence of the agitator is neglected and excess heat is removed by cooling water. The feed and product streams are present in the vapor phase. Molar balances for components A-H: 3

0 ) yi,6F6 - yi,7F7 +

∑ Vi,kRk k)1

(for i ) A, ..., H)

(30)

8092

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Energy balance for the reactor:

(

0 ) F6

)

H

3

∑yi,6Cpvi (T6 - Tr) - k)1 ∑ ∆Hr,kRk - Q˙ r i)A

(31)

0 ) yi,7F7 - yi,8(F8 + F9) - xi,10F10

H

∆Hr,k ) HoFk +



HiVi,k

numbers of moles Ni,r refer to the number of moles in the liquid phase only, because of the assumption that a buildup of these components in the vapor phase can be neglected. Product Separator. The molar balances for components A-H are given as

(32)

i)A

(for i ) A, ..., H) (45)

The energy balance for the separator is given as H

with

Hi )

yi,7Cpvi )(Tr - Ts) + HoVs - Q˙ s ∑ i)A

0 ) F7(

- Tref)

Cpvi (Tr

(46)

H

The reaction kinetics (eqs 5a-5c) remain the same for this model. Heat exchange with cooling water:

Q˙ r ) mcw,rCpcw(Tcw,r,out - Tcw,r,in) Q˙ r ) UAr

[

∆T1,r - ∆T2,r

ln(∆T1,r/∆T2,r)

]

HoVs )

v Ni,s RgTs pi,s ) Vv,s

(34) (35a)

∆T2,r ) Tr - Tcw,r,out

(35b)

(for i ) A, B, C)

sat (Ts) pi,s ) γi,sxi,10pi,s

(for i ) D, ..., H)

v Ni,r RgTr

(for i ) A, B, C)

Vv,r

(36) (37) xi,10 ) 0

)

(for i ) A, B, C)

xi,r )

∑ pi,r i)A

Ni,r

Ni,s H

(51)

(for i ) D, ..., H)

∑ Ni,s

(39) Vl,s )

(for i ) D, ..., H)

H

pi,s Ps

(52)

H

(for i ) A, ..., H)

(40)

(41)

∑ Ni,r i)D

l Ni,s ∑ i)D

(53a)

Fl,s

Vv,s ) Vs - Vl,s

(53b)

Heat exchange with cooling water:

Q˙ s ) mcw,sCpcw(Tcw,s,out - Tcw,s,in)

[

Q˙ s ) UAs

H



Vl,r )

(50)

i)D

H

pi,r yi,7 ) Pr

∑ pi,s

yi,8 ) yi,9 )

Bi sat pi,r (T) ) 10-3 exp Ai + Ci + Tr - Tref (for i ) D, ..., H) (38) Pr )

(49)

i)A

sat pi,r ) γi,rxi,rpi,r (Tr)

(

(for i ) D, ..., H)

(48)

H

Ps )

Vapor-liquid equilibrium:

pi,r )

(47)

Vapor-liquid equilibrium:

(33)

∆T1,r ) Tr - Tcw,r,in

∑xi,10F10Hvap,i

i)D

l Ni,r

i)D

(42a)

Fl,r

Vv,r ) Vr - Vl,r

(42b)

Reactor input stream (F6) and reactor output stream (F7):

F6 ) (0.8334 kmol s-1 MPa-1/2)xPm - Pr F7 ) (1.5355 kmol s

-1

MPa

-1/2

)xPr - Ps

(43) (44)

The constant values in eqs 43 and 44 are given by Jockenho¨vel et al.18 to match the base case. Note that the cooling water flux (mcw,r) is a direct control variable. Its dependence on a valve position is negligible. Note that, for components D-H, the

∆T1,s - ∆T2,s

ln(∆T1,s/∆T2,s)

]

(54) (55)

∆T1,s ) Ts - Tcw,s,in

(56a)

∆T2,s ) Ts - Tcw,s,out

(56b)

Compressor and Purge. The temperature changes that are due to the compressor work are taken into account by eq 57:

()

T8 ) Ts

Pm Ps

(1-κ/κ)

(57)

Stripper. This unit is modeled by a split fraction model, as proposed in Ricker and Lee.17 An energy balance was added by Jockenho¨vel et al.,18 and the split fractions (Φi) are modeled as third-degree polynomials in temperature. The pressures in the stripper and the mixing zone are assumed to be the same.

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8093 Table 2. Base Case Results

H

component

∑ (xi,10Cpli)(Ts - Tstr) + i)A

0 ) F10

Mole Fraction in Streams, yi,j 5

6

7

8

A B C D E F G H

0.45744 0.00448 0.43457 0.00124 0.07211 0.01197 0.02314 0.00630

0.31928 0.07872 0.24844 0.07024 0.19793 0.02431 0.04659 0.01449

0.27178 0.10765 0.18457 0.01030 0.18477 0.02702 0.13179 0.08212

0.31724 0.12249 0.22363 0.01508 0.20286 0.03372 0.06456 0.02042

flow (kmol/h) mass, SG/Ha

465.7 0.999104

1892.102

1476.0

1201.5

component

9

H

F4

(yi,4Cpvi )(T4 - Tstr) - HoVstr + Q˙ str ∑ i)A H

HoVstr )

∑ Hvap,i(yi,5F5 - yi,4F4)

flow (kmol/h) a

Q˙ str ) (2258.717 kJ/kg)m˘ cw,str

H

10

0.31724 0.12249 0.22363 0.01508 0.20286 0.03372 0.06456 0.02042

0.00000 0.00000 0.00000 0.00242 0.13427 0.02232 0.55138 0.28961

15.1

Vl,str )

11 0.00000 0.00000 0.00000 0.00024 0.00598 0.00104 0.54651 0.44623

259.5

Φi ) 1

211.3

yi,5 )

5

6

7

8

A B C D E F G H

0.43965 0.00453 0.46232 0.00072 0.05666 0.01939 0.02067 0.00826

0.32346 0.14777 0.18947 0.05933 0.16522 0.04123 0.04979 0.02373

0.28053 0.18403 0.11387 0.00587 0.15097 0.05199 0.1291 0.08363

0.32959 0.21655 0.13307 0.00798 0.15885 0.05467 0.0669 0.03239

flow (kmol/h) mass, SG/Ha

456.46 1.280171

2125.822

1707.081

1441.021

component

9

Mole Fraction in Streams, yi,j 0.32959 0.21655 0.13307 0.00798 0.15885 0.05467 0.0669 0.03239

flow (kmol/h)

9.40888

Ni,str (62)

Fstr

(for i ) A, B, C)

ai,j(Ts - 273)j ∑ j)0

(for i ) D, ..., H)

F5 ) F4 + F10 - F11

Mole Fraction in Streams, yi,j

A B C D E F G H

∑ i)D

10

11

0.00000 0.00000 0.00000 0.00140 0.10569 0.03622 0.48395 0.37273 253.563

0.00000 0.00000 0.00000 0.00014 0.00443 0.00160 0.60702 0.38682 210.885

(63)

3

Φi )

Table 3. Mode 1: 50/50 G/H Mass Ratio Results

a

(61)

Vapor-liquid equilibrium:

SG/H ) (yG,11MWG)/(yH,11MWH).

component

(60)

i)D

Mole Fraction in Streams, yi,j A B C D E F G H

(59)

xi,11 )

Φi(yi,4F4 + xi,10F10) F5

yi,4F4 + xi,10F10 - yi,5F5 F11 F

xi,11 ) (1 -

(64) (65)

(for i ) A, ..., H) (66)

(for i ) D, ..., F)

Ni,str

∑ xj,11) H j)D ∑ Nj,str j)D

(67)

(for i ) G, H) (68)

Step 2.7: Model Solution. Because of the complexity of the system, the information generated from stage 1 will be helpful to guide/identify the solution’s location of the system. With the modifications made by Jockenho¨vel et al.,18 there are 30 states, and 11 manipulated variables, compared to the model in Ricker and Lee,19,17 which has 26 states and 10 manipulated variables. As it was done in steps 1.6 and 1.7, the model equations (eqs 26-68) are presented at steady state in a vectorial form; thus,

SG/H ) (yG,11MWG)/(yH,11MWH).

0 ) f[x,u,p]

(69)

x ) [xs,xp]

(70)

The heating medium is saturated steam, which condenses completely at a constant temperature. The enthalpy of the steam has been chosen to fit the steam flux, with the heat duty of the stripper given by Downs and Vogel15 for the base case. Molar balances for components G and H:

xs ) [Tm, Tr, Ts, Tstr, Ni,m, Ni,r, Ni,s, NG,str, NH,str]T (for i ) A, ..., H) (70a)

0 ) (1 - Φi)(xi,10F10 + yi,4F4) - xi,11F11

xp ) [F5, F6, F7, Pm, Pr, Ps, T8, Tcw,r,out, Tcw,s,out, Vl,r, Vv,r, Vl,s, Vv,s, Vl,str, Hi, Ho Vs, HoVstr, Q˙ r, Q˙ s, Q˙ str,∆T1,r, ∆T2,r, (for i ) G, H) (58)

Energy balance for the stripper:

where

sat sat , pi,s, pi,s , xi,r, xi,10, xi,11, yi,5, yi,6, yi,7, ∆T1,s, ∆T2,s, pi,r, pi,r

yi,8, yi,9, Rk, ∆Hr,k, Φi,]T (for i ) A, ..., H; k ) 1, ..., NR) (70b)

8094

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

u ) [F1, F2, F3, F4, F8, F9, F10, F11, mcw,r, mcw,s, mcw,str]T

(71)

p ) [T1, T2, T3, T4, Tcw,r,in, Tcw,s,in, Vm, Vr, Vs, UAr, UAs, Rg, κ, Fr, Fs, Fstr, Ai, Bi, Ci, Cpli, Cpvi , Cpcw, HoFk, Hvap,i, xi,1, (for c ) D, ..., H; xi,2, xi,3, xi,4, γi,r, γi,s, k, ac,j]T i ) A, ..., H; j ) 0, ..., 3; k ) 1, ..., NR) (72) The state vector x is comprised of 30 state variables (xs) and 129 process variables (xp); vector u contains 11 design variables, and 100 properties and parameters of the system are listed in vector p (Tables A1-A5 in the Appendix provides the values for these variables). The system consists of 270 variables and 159 equations. Thus, NDOF ) 111. Therefore, by specifying vector u and p for different process specifications, the system can be solved. A sequential modular approach has been used to solve the steady-state model (see eq 69). Here, the selection of dependent and independent variables becomes simpler,20 if it is assumed that all equipment models are written as simulation models, i.e., the output flows are calculated from given input flows, conditions, and equipment dimensions. All input flows and equipment dimensions are independent variables, specified by the “user”, and all equipment outlet streams are dependent variables. However, when recycle flows appear, it is necessary to apply tearing techniques21 to initialize the flowsheet appropriately. The common selection when only one loop is encountered is the stream entering the reactor, which, in this case, is stream 6. Thus, the corresponding module (unit) calculation sequence is

Figure 12. G/H mass ratio, as a function of feed inlet ratio (fC,4).

Reactor f Separator f Stripper f Mixing Zone After solution of each unit, the connectivity information (i.e., the inlet flow rate, pressure P, and temperature T) becomes available for the following unit. Step 2.8: Model Analysis. Step 2.8.1. Analysis. Three different operation modes were analyzed in this process. Table 1 lists the manipulated variables (vector u) corresponding to these modes.17 In this respect, the results obtained from the modular approach, in terms of the compositions and flow rates of the streams, are given in Tables 2-4. On the other hand, the energy balances have been included for the purpose of being able to use a model that is able to reproduce the behavior of the process in a more-realistic manner. Thereby, Table A6 in the Appendix shows the solutions of the state vector xs obtained with the non-isothermal model (this work, based on the reported research by Jockenho¨vel et al.18) and the reported values for the isothermal model22 for the three operation modes analyzed. As can be observed, some of the major differences are in the product separator, where the vaporliquid equilibrium has a key role, as well as the amount of liquid that is vaporized. These differences can be assigned mainly because of the addition of energy balances to the process model. Step 2.8.2. Identification of Operational Targets. Steps 2.1 and 2.2 identify the operating targets. However, for this problem, the design targets are defined in advance (see Table 1 in Downs and Vogel15), i.e., a G/H mass ratio of 50/50 (the base-case scenario). Step 2.9: Simulations Agreement. Figure 12 compares the results obtained with the simple (stage 1) and detailed (stage 2) models, as a function of the ratio of component C, with respect to component A in stream 4 (fC,4 (yC,4/yA,4)). As can be observed, both curves follow the same trend, although the mass

Figure 13. Reactor holdup profiles for the base case. Table 4. Mode 2: 10/90 G/H Mass Ratio Results Mole Fraction in Streams, yi,j component

5

6

7

8

A B C D E F G H

0.42958 0.00443 0.45172 0.00005 0.08549 0.02776 0.00399 0.01527

0.34465 0.07869 0.19418 0.01153 0.26240 0.05793 0.00920 0.04142

0.30864 0.09704 0.12464 0.00107 0.22146 0.07249 0.02492 0.14974

0.36252 0.11352 0.14587 0.00059 0.23337 0.07610 0.01222 0.05581

flow (kmol/h) mass, SG/Ha

430.4475 0.110807

2062.893

1683.037

1413.116

component

9

10

11

A B C D E F G H flow (kmol/h)

0.36252 0.11352 0.14587 0.00059 0.23337 0.07610 0.01222 0.05581 16.09766

0.00000 0.00000 0.00000 0.00010 0.15672 0.05097 0.09458 0.69763 243.825

0.00000 0.00000 0.00000 0.00001 0.00726 0.00246 0.11842 0.87184 194.638

Mole Fraction in Streams, yi,j

a

SG/H ) (yG,11MWG)/(yH,11MWH).

selectivity values (measured using the mass ratio SG/H) obtained from the two models do not match. This is not surprising, because these differences are due to the simplifications made in the stage 1 model, particularly, in the separation section, where the effects of the vapor/liquid separator and stripper are

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8095

Steps 3.3. and 3.4: Simulation and Verification of the Results. This particular case study exhibits unstable performance in the open loop. Consequently, it is important to analyze the dynamic behavior of the system. In this respect, the nonlinear model represented by eq 69 is presented now in its dynamic form:

dxs ) x3 s ) f[x,u,p] dt

(73)

y ) h[x,u,p]

(74)

where the vectors x, u, and p remain unchanged in eqs 70, 71, and 72, respectively, and vector y represents the output variables (variables to be measured), which is a function h of the state variables x, the manipulated variables vector u, and system parameters vector p. Equation 75 lists the variables considered in the output vector y. Figure 14. Temperature and pressure profiles for the base case.

lumped into a “simple” separation factor (Ri,7); thus, the thermodynamics effects and inherent nonlinearities were hidden. Nevertheless, the stage 1 model is considered, to some extent, to be representative of the behavior of the process; therefore, the decision was made to proceed further with the analysis of the detailed model. The flow ratio fC,4 is considered to be a process disturbance (IDV(1)15); therefore, it is important to know its effect on the performance of the system. Point “A” in Figure 12 is the nominal operating point (the base-case scenario) for the aforementioned flow ratio. Figure 12 shows that, if there was the interest to increase the selectivity of G to H, there would exist two possibilities at two different fC,4 values, although two different production rates will also be achieved. Step 2.10: Storage of Information. At this point of the model-based Methodology, it is possible to determine the operating conditions for the base-case scenario from the stage 2 model. The sensitivity analysis from the stage 1 models revealed patterns, in terms of manipulated and process variables, that were verified using more-detailed models (see Figure 12). Thus, the information generated has helped us to recognize its behavior and characterize it in terms of the flow rate, P, and T values in all the streams. 3.4. Stage 3: Verification. The objective of this stage is that, from the development and analysis of simple and detailed models, it would become easier to validate the actual process behavior through the use of more-rigorous models, either by dynamic or steady-state simulations. Another objective is that, from the information generated by the stage 1 and 2 models, it would be possible to outline an appropriate control structure for a given design target(s), because of the static responses of the system are known from the analysis of sensitivity performed. Step 3.1: Flowsheet Setup. Figure 1 remains unchanged. For the TE problem, for proprietary reasons, the amount of information is limited, making the use of a commercial process simulator infeasible. Step 3.2: Generated Knowledge from Previous Stages. Because of the aforementioned reasons, the model implemented in the stage 2 model will be considered to be the “rigorous” model, because it reproduces the data from the literature. On the other hand, from the analysis in the stage 1 and 2 models, the design decisions have also been verified using static models.

y ) [F11,Pr,Tr,Tstr,Vl,r,Vl,s,Vl,str,yB,9,yC,9]T

(75)

Consequently, the dynamic model has 30 ordinary differential equations (ODEs) (states) and 129 algebraic equations (AEs), with 11 control variables. Figures 13 and 14 show the (open-loop) dynamic response for the reactor holdups, as well as the temperature, liquid volume, and pressure profiles in the equipment, respectively. Small oscillations are observed with increasing amplitude in the profiles during the first hour of operation. After this initial period, the amplitude of the oscillations increases substantially, mainly in the temperature of the reactor. Hence, the instability of the process is confirmed, because the safety limits have been reached and, therefore, the operation must be shut down. In Figure 15, the phase planes of the component A and G holdups, with respect to the reactor temperature, are presented. The limit cycles that are observed are a clear indication of the (open-loop) unstable nature of the process. Within the first hour of operation, the system oscillates around a point relatively close to the steady-state value, but, with increasing time, the number of cycles increases considerably. Now, if eqs 73 and 74 are linearized around a given operating point (ss), then

∂f dx ∂f ) | (x - xss) + |ss(u - uss) dt ∂x ss ∂u y)

∂h | (x - xss) ∂x ss

Therefore,

dx ) A(x - xss) + B(u - uss) dt y ) C(x - xss) where A is an n × n state matrix, B is an n × r input matrix, and C is an m × n output matrix, all of which have been evaluated at a given steady state for n states, r inputs, and m outputs. Therefore, another test to verify the stability of a process is the analysis of the eigenvalues (λi) of the Jacobian of the plant (matrix A) evaluated at steady state. In Table 5, the eigenvalues for matrix A for the base-case scenario are listed. In this case, five positive values appear: the pair-conjugate λ14λ15, λ23, λ24, and λ26. These results are an indication of the

8096

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Figure 15. Limit cycles of the base-case scenario: (a) NA,r vs Tr and (b) NG,r vs Tr.

Figure 16. Control structure proposed for the TE problem. Table 5. Eigenvalues of A for the Base Case parameter

value/expression

parameter

value/expression

parameter

value/expression

λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10

-0.42525 -0.20689 -0.09266 -0.04290 -0.03216 + 0.00011i -0.03216 - 0.00011i -0.02236 -0.01941 + 0.00671i -0.01941 - 0.00671i -0.01625 + 0.00169i

λ11 λ12 λ13 λ14 λ15 λ16 λ17 λ18 λ19 λ20

-0.01625 - 0.001685i -0.00807 -0.00806 0.00106 + 0.00369i 0.00106 - 0.00369i -0.00396 + 0.00128i -0.00396 - 0.00128i -0.00404 -0.00371 -0.00299

λ21 λ22 λ23 λ24 λ25 λ26 λ27 λ28 λ29 λ30

-0.00176 -0.00065 0.00044 0.00041 -5.382 × 10-5 6.125 × 10-6 -3.753 × 10-6 -3.845 × 10-5 -2.227 × 10-5 -2.164 × 10-5

Table 6. Control Structure Proposal for TE Problema uS1

a

meaning

Da mcpcw Ri,7

reaction rate, wrt feed rate (process) reactor heat removal recovery of i from separation section

fi σ

feed flow ratio purge factor

k(f(Pr)), wrt FA,1 mcw,r, wrt FA,1 Fi,rec, wrt Fi,r Fi,in, wrt FA,1 F9, wrt Ftop,s

S1 denotes stage 1; S2 denotes stage 2. b Data taken from Luyben et al.7

instability of the process, which also has been shown using the dynamic simulations. With respect to the operability of the system, it is necessary that all unstable states are observable (i.e., the plant is detectable) and that these states also are controllable (i.e., the plant is stabilizable). However, such analysis has not been conducted within this work. Step 3.5: Storage of Information. Several tests have been performed (open-loop simulations, eigenvalues analyses, con-

c

uS2

yS2

yb

mcw,r mcw,r mcw,s mcw,str FD,2/FE,3 F9 F11 F10 F8 FC,4

Pr Tr Vl,s Tstr Vl,r yB,9 F11 Vl,str yields yC,9

yA,9 Tr Vl,sc Tstr Vl,r yB,9 F11 Vl,str yields Pr

Luyben et al.7 used the cooling water flow to the condenser.

trollability matrix analyses) to validate the unstable nature of the TE problem. It also has been verified that, as noted by Downs and Vogel15 and Ricker and Lee,17 after the first hour of operation, the process starts to become unstable, reaching the safety limits; the operation then must be shut down at that point. From the previous analysis, it is possible now to relate the different effects that are encountered, to outline a control

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8097

Figure 17. Dynamic response for output variables in a closed loop, based on the control structure proposal for the base case.

structure that leads to a feasible control strategy. This is done in the following discussion. 3.4.1. Control Structure Implications. The purpose of this section is to try to relate the insights obtained from the models developed in the stage 1 and 2 models and their corresponding analysis with a sketch of a control structure. While the following description relies on the procedure developed by Luyben et al.,7 for the synthesis of a control structure, it will become more clear that the selection of the manipulated and controlled variables can be facilitated through the information (insights) generated using the model-based Methodology. It is also intended to show how the DOFs in the control problem can be fulfilled with the knowledge obtained through the model-based Methodology. The plant Da value (eq 6) assesses the ratio of the rate of reaction, with respect to the feed to the process of one of the limiting reactants. In the analysis that has been conducted, in the calculation of Da, the parameters Pr and Vr were assumed to be kept constant, as well as the temperature at the inlet of the reactor (Tm). Therefore, a variation in the Da value can be interpreted, mainly, as a change in the feed of A to the plant. Ricker and Lee17 assumed that Tm remained constant through the calculations and scenarios analyzed. In this context, for this dimensionless group, the manipulated variables would be FA,1 and the variables to be controlled would be Pr. The variable mcpcw relates the feed of cooling water to the reactor, with respect to the feed of component A to the process. The analysis in this dimensionless variable reveals that the heat removal of the reactor (Q˙ r) is affected by changes in Tcw,r,out. Lyman and Georgakis23 proposed to control Tcw,r,out via mcw,r in an inner loop and, the master loop Tr would be controlled by the manipulation of Tcw,r,out,sp. It is logical to expect that the more heat that is removed from the reactor, the lower the selectivity that is obtained.

In stage 1, the separation was modeled through the recovery factors Ri,7, which helps to determine the process conditions to achieve a given set of specifications. Note that the recovery factors are defined as the ratio of flow rate of component i in the recycle stream with respect to the flow rate of component i from the reactor. Therefore, their control implications can be understood as controlling the inventories through the plant to satisfy the product specifications. For example, in the control structure proposed by Luyben et al.,7 the gas recycle (F8) is used to control the yields; Lyman and Georgakis23 used a temperature controller to manipulate the steam flow rate to control the composition of E at the bottom of the stripper. Both Luyben and Lyman used the purge flow rate (F9) to maintain the composition of B (the former in the purge stream and the latter in the recycle stream). On the other hand, because the recycle stream is in the gas phase, the possibility to control either the pressure or temperature through the cooling water in the separator or through the steam in the stripper is considered. The analysis of the flow ratio helps to evaluate the process response under different feeds of reactants (disturbance) to the system. It can also help to determine appropriate setpoints needed to maintain a given operating condition, such as composition or flow rate. In step 6 of Luyben’s procedure, one is required to fix a flow in every recycle loop and then select the best manipulated variables to control inventories (pressures and levels). Based on this information, Luyben et al.7 selected the ratio FD,2/FE,3 ((yD,2F2)/(yE,3F3)) to control the liquid level in the reactor, as well as to maintain the desired product distribution G and H. In the analysis performed in this case study, the flow ratio considered was mainly with respect to component A. This analysis revealed that, in the case of fC, as the ratio was increased, the selectivity also increased. This is not surprising, because the greater the amount C fed, the larger the production

8098

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

Table A1. Physical Properties of the Componentsa

Table A3. Kinetic Values in Rate Expressions

Heat Capacity (kJ kg-1 °C-1) molecular liquid density, component weight, MW F (kg/m3) liquid A B C D E F G H

2.0 25.4 28.0 32.0 46.0 48.0 62.0 76.0

299.0 365.0 328.0 612.0 617.0

vapor

heat of vaporization (kJ/kg)

14.6 2.04 1.05 1.85 1.87 2.02 0.712 0.628

7.66 4.17 4.45 2.55 2.45

parameter

value

parameter

value

parameter

value

c1,1 c2,1

44.06 42600

c1,2 c2,2

10.27 19500

c1,3 c2,3

59.50 59500

Table A4. Dimensionless Variables 202 372 372 523 486

Vapor Pressure Constants, Using the Antoine Equationb component

A

B

C

D E F G H

20.81 21.24 21.24 21.32 -3318.0

-1444.0 -2114.0 -2144.0 -2748.0 250

259 266 266 233 486

parameter

value

γ1 γ2 γ3 δc T2 K* 2 K* 3 mcpcw

59.67809 27.37744 83.35321 74.45931 359.25 K 1.763731 3.11015 × 10-5 62.62184

Table A5. Parameter Specifications for TE Problema

a Data taken from ref 15. The Antoine equation is given as P ) exp[A i + Bi/(Ci + T(°C))]. Parameter values for components A, B, and C are not listed, because they are effectively noncondensible.

Table A2. Stream Conditions for the Base Casea Mole Fraction in Streams, yi,j component A B C D E F G H flow temp, T a

1

2

3

4

0.99990 0.00010 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

0.00000 0.00010 0.00000 0.99990 0.00000 0.00000 0.00000 0.00000

0.00000 0.00000 0.00000 0.00000 0.99990 0.00010 0.00000 0.00000

0.48500 0.00500 0.51000 0.00000 0.00000 0.00000 0.00000 0.00000

11.2 kmol/h 45.0 °C

114.5 kmol/h 45.0 °C

98.0 kmol/h 45.0 °C

417.5 kmol/h 45.0 °C

Data taken from ref 15.

of G. Consequently, the heat duty increases, because the first reaction is the most exothermic. In this respect, the heat duty in the reactor can be considered to be an indirect indication of changes in the level of the reactor. Given the aforementioned considerations, Table 6 shows the control structure proposed for the TE problem. This table shows that most of the DOFs for the control problem are covered through the model-based analysis of the process. However, the remaining DOFs to be fulfilled follow a logical sequence. For instance, the flow rate product (F11) is actually the control objective: then this is set on flow control. The control of the level of the stripper (Vl,str) can be paired with the bottom flow from the separator (F10), because the steam flow to the stripper has already been used to control the temperature (Tstr). Douglas24 recommended that the valve be kept wide open (maximizing the flow) in gas recycle processes, to improve the yields. Therefore, one DOF is removed. The remaining variable (FC,4, given that it is the largest gas feed) will be used to control the inventory of C through the system. Figure 16 shows the control structure applied to this system. Based on this control structure proposal (see Table 6), Figure 17 shows the dynamic trajectories in the closed loop for the (controlled) output variables (y). The simulations were performed only for the base-case scenario. It is important to remark that only proportional controllers have been used to control the

a

Data taken from ref 18.

system. As can be seen, the output variables (y) result adequately controlled based on this implementation, thereby, verifying the insights generated from the model-based Methodology. 4. Conclusions The present case study has demonstrated how, through the model-based analysis, it is possible to start from a conceptual design to the point where the implementation of a control strategy is devised to achieve (or maintain) a given set of design targets. The Tennessee-Eastman (TE) problem represents a very interesting problem of integration, in that, through the model-based Methodology, it is possible to identify some of the underlying phenomena involved. Most of the research conducted for the TE problem mainly addresses the control problem. However, little attention has been paid to designrelated issues, which are a necessary step to perform successive control studies. On the other hand, although the control

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8099 Table A6. Comparison of the Results between the Modular Approach and that of Ricker and Lee17 Base Case holdup

this work

NA,m NB,m NC,m ND,m NE,m NF,m NG,m NH,m Tm NA,r NB,r NC,r ND,r NE,r NF,r NG,r NH,r Tr NA,s NB,s NC,s ND,s NE,s NF,s NG,s NH,s Ts NG,str NH,str Tstr

47.236 11.646 36.755 10.392 29.283 3.097 6.093 2.143 368.367 4.296 1.702 2.918 0.191 11.733 1.520 76.388 71.486 394.791 26.333 10.167 18.563 0.069 4.173 0.627 15.495 11.138 355.704 21.492 17.549 338.859

Mode 1 (50/50)

Ricker and

Lee17

48.838 13.493 40.030 10.442 28.488 2.514 5.403 2.517 359.250 3.972 1.665 2.888 0.169 10.641 1.295 67.260 68.196 395.166 28.895 12.119 21.022 0.097 5.905 0.719 20.492 16.195 353.250 21.754 17.510 337.863

this work 54.379 24.767 31.779 9.759 27.465 6.634 8.823 4.361 388.951 4.593 3.006 1.860 0.112 9.857 3.353 76.667 75.049 398.272 28.974 18.986 11.634 0.059 4.275 1.693 22.852 18.032 368.281 27.766 17.698 339.621

structure proposed turned out to be similar to Luyben et al.,7 the synthesis of the control structure was intended to follow a more-fundamental basis, because it was performed through the model-based analysis, instead of from a priori knowledge. With this example, it was also demonstrated how the modelbased analysis helps to verify the design decisions that must be made in the modeling stages. In addition, it also demonstrated how the design and control problems can be related (integrated) simultaneously, which is an ultimate objective of this work. 5. Appendix Table A1 lists the physical properties for the components in the Tennessee-Eastman (TE) problem. Table A2 lists the stream conditions utilized for the basecase scenario.15 The molar compositions and temperature do not change between operation modes. Streams 1-4 are combined to produce stream 1 in the simplified flowsheet (see Figure 2). Table A3 gives the parameter values used in eqs 5a-5c. Table A4 gives the values for the dimensionless variables used in the stage 1 model when the energy balance is added to the model, based on the base-case scenario and the definitions given in step 1.7. Table A5 gives the parameter specifications (vector p), and Table A6 shows the solutions of the state vector xs obtained with the non-isothermal model (from this work, based on Jockenho¨vel et al.18) and the reported values for the isothermal model22 for the three operation modes analyzed. Nomenclature Cpi ) heat capacity of component i (kJ kmol-1 K-1) Fj ) molar flow rate of stream j (kmol/s)

Ricker and

Mode 2 (10/90) Lee17

52.847 24.189 30.924 9.729 26.940 6.747 8.125 3.870 370.523 4.315 2.834 1.750 0.106 9.025 3.117 68.626 65.786 398.357 31.005 20.415 12.487 0.080 5.967 2.057 28.442 22.117 368.971 27.757 17.716 339.635

this work

Ricker and Lee17

56.721 12.951 31.958 1.898 43.185 9.534 1.515 6.817 366.052 7.299 2.409 2.951 0.012 8.233 2.763 8.539 80.491 393.684 36.724 11.500 14.777 0.007 4.770 1.576 3.284 23.880 367.284 4.351 32.035 338.456

59.874 14.107 33.449 2.206 43.729 9.645 1.708 7.440 359.250 6.090 1.956 2.432 0.016 9.972 3.285 10.927 95.896 397.350 32.039 10.291 12.794 0.008 5.680 1.871 4.032 28.679 363.450 4.363 32.037 338.550

Hi ) enthalpy of component i (kJ) HoFk ) enthalpy of formation of reaction k (kJ/kmol) HoVj ) enthalpy of vaporization in equipment j (kW) Ni,j ) holdup of component i in equipment j (kmol) Da ) Damko¨hler number Pj ) total pressure in equipment j (kPa) pi,j ) partial pressure of component i in equipment j (MPa) Q˙ j ) heat duty in equipment j (kW) Rg ) gas constant (kJ kmol-1 K-1) Rk ) reaction rate of reaction k (kmol h-1 m-3) Tj ) temperature of stream j (K) UAj ) specific heat-transfer rate coefficient of equipment j (kW/ K) Vj ) volume of equipment j (m3) xi,j ) mole fraction in the liquid of component i in stream j Yi ) yield of component i yi,j ) mole fraction in the vapor of component i in stream j Greek Symbols Ri,j ) recovery factor of component i, with respect to stream j γi,j ) activity coefficient of component i in stream j γk ) Arrhenius number in reaction k k ) reaction constant (kmol h-1 m-3) κ ) compression factor λ ) purge factor νi,k ) stoichiometric coefficient of component i in reaction k σ ) purge relationship; σ ) 1/(1 - λ) Φi ) split factor for component i in the stripper Subscripts cw ) cooling water i ) component index j ) stream index k ) reaction index l ) liquid phase

8100

Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007

m ) mixing zone r ) reactor s ) separator str ) stripper v ) vapor phase Superscripts l ) liquid phase sat ) saturation conditions v ) vapor phase Acknowledgment The authors gratefully acknowledge the financial support from El Consejo Nacional de Ciencia y Tecnologı´a (Me´xico) for the realization of this study. Literature Cited (1) Ramı´rez, E.; Gani, R. Methodology for the Design and Analysis of Reaction-Separation Systems with Recycle. 1. The Design Perspective. Ind. Eng. Chem. Res. 2007, 46, 8066-8083. (2) Larsson, T.; Skogestad, S. Plantwide controlsA review and a new procedure. Model., Identif. Control 2000, 21, 209-240. (3) Nishida, N.; Stephanopoulos, G.; Westerberg, A. W. A Review of Process Synthesis. AIChE J. 1981, 27, 321-351. (4) Fisher, W. R.; Doherty, M. F.; Douglas, J. M. The Interface between Design and Control. 1. Process Controllability. Ind. Eng. Chem. Res. 1988, 27, 597-605. (5) Fisher, W. R.; Doherty, M. F.; Douglas, J. M. The Interface between Design and Control. 2. Process Operability. Ind. Eng. Chem. Res. 1988, 27, 606-611. (6) Fisher, W. R.; Doherty, M. F.; Douglas, J. M. The Interface between Design and Control. 3. Selecting a Set of Controlled Variables. Ind. Eng. Chem. Res. 1988, 27, 611-615. (7) Luyben, W. L.; Tyre´us, B. D.; Luyben, M. L. Plantwide Process Control; McGraw-Hill: New York, 1998. (8) Skogestad, S. Control structure design for complete chemical plants. Comput. Chem. Eng. 2004, 28, 219-234. (9) Kiss, A. A.; Bildea, C. S.; Dimian, A. C.; Iedema, P. D. Design of Recycle Systems with Paralell and Consecutive Reactions by Nonlinear Analysis. Ind. Eng. Chem. Res. 2005, 44, 576-587. (10) Luyben, M. L.; Floudas, C. A. Analyzing the Interaction of Design and Controls1. A Multiobjective Framerwork and Application to Binary Distillation Synthesis. Comput. Chem. Eng. 1994, 18, 933-969.

(11) Luyben, M. L.; Floudas, C. A. Analyzing the Interaction of Design and Controls2. Reactor-Separator-Recycle System. Comput. Chem. Eng. 1994, 18, 971-994. (12) Mohideen, M.; Perkins, J.; Pistikopoulos, E. A Framework for Process Design and Control. In Control ’96, UKACC International Conference on (Conf. Publ. No. 427), Vol. 2; Institution of Electrical Engineers (IEE): London, U.K., 1996. (13) Mo¨nnigmann, M.; Marquardt, W. Parametrically robust controlintegrated design of nonlinear systems. In American Control Conference, 2002. Proceedings of the 2002, Vol. 6; American Automatic Control Council: Anchorage, AK, 2002. (14) Downs, J.; Vogel, E. A plant-wide industrial process control problem. Presented at the AIChE Annual Meeting, Chicago, IL, 1990. (15) Downs, J.; Vogel, E. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245-255. (16) Glasser, D.; Hildebrandt, D.; Crowe, C.; A Geometric Approach to Steady Flow Reactors: The Attainable Region and Optimization in Concentration Space. Ind. Eng. Chem. Res. 1987, 26, 1803-1810. (17) Ricker, N.; Lee, J. Nonlinear modeling and state estimation for the Tennessee Eastman Challenge Process. Comput. Chem. Eng. 1995, 19, 983-1005. (18) Jockenho¨vel, T.; Biegler, L. T.; Wa¨chter, A. Dynamic optimization of the Tennessee Eastman process using the OptControlCentre. Comput. Chem. Eng. 2003, 27, 1513-1531. (19) Ricker, N.; Lee, J. Nonlinear model predictive control of the Tennessee Eastman Challenge Process. Comput. Chem. Eng. 1995, 19, 961981. (20) Wells, G.; Rose, L. The Art of Chemical Process Design; ComputerAided Chemical Engineering, Vol. 2; Elsevier: Amsterdam, 1986. (21) Sargent, R.; Westerberg, A. “Speed-Up” in Chemical Engineering Design. Trans. Inst. Chem. Eng. 1964, 42, T190-T197. (22) Ricker, N. Optimal Steady-State Operation of the Tennessee Eastman Challenge Process. Comput. Chem. Eng. 1995, 19, 949-959. (23) Lyman, P.; Georgakis, C. Plant-wide control of the Tennessee Eastman Problem. Comput. Chem. Eng. 1995, 19, 321-331. (24) Douglas, J. M. Conceptual Design of Chemical Processes, 1st Edition; International Series of Chemical Engineering; McGraw-Hill: New York, 1988.

ReceiVed for reView October 16, 2006 ReVised manuscript receiVed August 3, 2007 Accepted August 8, 2007 IE061328P