An Optimization-Based Approach for the Operability Analysis of

Methodology for the Steady-State Operability Analysis of Plantwide Systems. Sivakumar ... Optimal operability by design in a methanol reforming-PEM fu...
4 downloads 4 Views 171KB Size
4238

Ind. Eng. Chem. Res. 2001, 40, 4238-4252

An Optimization-Based Approach for the Operability Analysis of Continuously Stirred Tank Reactors Sivakumar Subramanian, Derya Uztu 1 rk, and Christos Georgakis* Chemical Process Modeling and Control Research Center and Department of Chemical Engineering, Lehigh University, Bethlehem, Pennsylvania 18015

A generalized methodology for examining the inherent operability characteristics of chemical reactors is presented. This approach is developed by extending the operability framework proposed by Vinson and Georgakis (In DYCOPS-5, 5th IFAC Symposium on Dynamics and Control of Process Systems; Pergamon Press, 1998 and J. Process Control 2000, 10, 185-194) and Uztu¨rk and Georgakis (Ind. Eng. Chem. Res. 2001, manuscript submitted) to nonsquare systems. This extension is motivated by the fact that operability problems with more control variables than set-point-controlled outputs are more common. In this formulation, both the steady-state and dynamic operability analyses require the solution of families of optimization problems. The proposed methodology is employed to investigate the integration of the design and control of systems of a single CSTR and two CSTRs in series. Results obtained show that (i) reactors operating at higher temperatures are inherently faster than reactors operating at lower temperatures, (ii) systems with two reactors in series are faster than a single reactor operating at the same temperature, and (iii) two reactors in series designed to have the same heat removal load respond much faster than all other systems studied. 1. Introduction Operability analysis provides an opportunity for process designers to address control issues quite early in the design process. The importance of such evaluations has been increasingly emphasized in the literature in the past decade or so. Ever-increasing competition in the global market and increasing environmental regulations force chemical processing industries to operate with the maximum efficiency and the minimum emissions possible. These factors often lead to process designs with considerable process integration. It is wellknown that such designs can pose difficult control problems. In such a scenario, it is imperative to test every prospective plant design for its operability. To effectively accomplish this task, the designer needs a metric that can be used to measure the operability of a given design. As discussed by Morari and Perkins1 and Vinson and Georgakis,2,3 such a measure has been a missing link in the integration of process design and control. Some of the earlier attempts to bridge this gap for a nonlinear plant include: flexibility analysis,4,5 backoff from optimum,6,7 and bifurcation analysis.8 In flexibility analysis, an index is defined as the size of the largest normalized hypercube that fits in the uncertainty region without violating the process constraints. The backoff approach, on the other hand, prepares the designer for the worst disturbance by moving the optimum operating point a safe distance away from the constraints so that the plant would be operable even if disturbances entered the plant. The bifurcation-based approach delineates the regions of multiplicity in the parameter space and, thus, cautions the designer to avoid them if possible. * Corresponding author. Address: CPMC Research Center, D-311 Iacocca Hall, 111 Research Drive, Lehigh University, Bethlehem, PA 18015. Phone: (610) 758-5432. Fax: (610) 7585297. Email: [email protected].

Recently, Vinson and Georgakis2,3 proposed a simple yet powerful metric to quantify the operability of processing systems. The technique has proven to be effective for both linear and nonlinear9,10 processes. Its extension to dynamic operability analyses is being studied by Uztu¨rk and Georgakis.11 In these works, for ease of development and exposition, only process models with equal numbers of inputs and outputs were considered. However, many realistic problem definitions would have more inputs than equality specifications on the outputs. In general, on the basis of operational requirements, process outputs can be classified into two broad categories: (1) set-point controlled, outputs to be controlled at a desired value, and (2) set-interval controlled, outputs to be controlled within a desired range. For instance, production rate and product quality might fall into the first category, whereas process variables, such as level, pressure, and temperature in different streams/ units, might fall into the second category. One of the objectives of the present work is to extend this operability framework to problems with extra degrees of freedom, that is, processes with more input variables than set-point-controlled outputs. This article demonstrates the proposed approach for the case of designing CSTRs in which a first-order, irreversible, exothermic reaction of the type A f B is taking place. As Russo and Bequette8 remark, such CSTR systems are, perhaps, the most extensively studied in chemical engineering literature. Surprisingly, not many articles have addressed the issue of interaction between their design and control characteristics, motivating the present work. The designs studied in this work include a single CSTR and a series of two CSTRs. Steady-state design procedures for a single reactor are standard. In the design of a system with two reactors in series, several alternatives arise because of the extra degrees of freedom in the design. We explored many alternatives by using two reactors operating at different

10.1021/ie001111+ CCC: $20.00 © 2001 American Chemical Society Published on Web 06/08/2001

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4239

temperatures, as well as by using reactors with different volume ratios. A selected set of these CSTR systems was subjected to rigorous steady-state and dynamic operability analyses. Luyben12 studied the CSTR design problem with an emphasis on operability issues. He used the Qmax/Q ratio as a measure of controllability, where Q is the nominal heat load of a design and Qmax is the heat that can be removed with the maximum coolant flow rate when the reactor is operating at its nominal temperature. He argues that a higher Qmax/Q ratio implies better controllability. When single-CSTR designs operating at different temperatures are compared, he found that the Qmax/Q ratio can go through a maximum, leading to two design temperatures with the same Qmax/Q ratio. In that case, one would obviously opt for the reactor design operating at higher temperature because of economic factors. He also studied the closed-loop responses of single-reactor and reactors-in-series configurations using PI controllers. He found that a single CSTR copes with large disturbances very well, whereas multiple CSTRs in series produce large overshoots. Jaisathaporn13 analyzed the load rejection and switchability (ability to go from one set point to another) characteristics of three single-CSTR designs. The time spent in producing off-specification products was used as a measure of lack of controllability. He found that large reactors have better load rejection capabilities. On the other hand, smaller reactors showed fast dynamics during set-point transitions, but their responses were oscillatory, producing more off-specification products. Therefore, on both counts, he recommended using larger reactors. In studies of this nature, controllability characteristics are intricately linked to the controller choice, configuration, pairing, and tuning employed. Consequently, they do not reveal the inherent operability characteristics of the process. The minimum-time formulation adopted from Uztu¨rk and Georgakis11 for this work evaluates the operability of the process without considering a preassigned feedback controller. Although this is an open-loop controller, it establishes a theoretical upper bound on the performance of any feedback controller and identifies the inherent limitations of the process. Douglas and Denn14 and Douglas15 presented time-optimal control studies for a single CSTR with the purpose of developing multivariable controllers. They derived an analytical solution for this optimal control problem considering a simple two-state model of the CSTR process. In this work, we utilized a more realistic model of the process for which the time-optimal control problem is formulated as a nonlinear programming problem and solved numerically as described in Appendix B. Although we study this rather simple reacting system, the methodology presented is equally applicable to more complex problems as well. However, such applications would require the solution of larger optimization problems. The outline of the paper follows. In the next section, we briefly discuss the operability index of Vinson and Georgakis2,3 and extend this methodology to nonsquare systems. Then, in section 3, we present the process model and the specifications of the CSTR process. In sections 4 and 5, we present results of the steady-state and dynamic operability analyses of single reactors and systems with two reactors in series, respectively. Section 6 presents the important conclusions of the paper. Appendices A and B document the detailed process

model and model parameters and the details of the optimal control formulation, respectively. 2. Process Operability Analysis In this section, we first review the steady-state and the dynamic operability techniques, which are based on the operability analysis framework introduced by Vinson and Georgakis.2,3 Next, some extensions and algorthmic developments are described. In this report, we refer to the following definition of operability adopted from Vinson:16

Operability. A process is operable if the available set of inputs is capable of satisfying the desired steady-state and dynamic performance requirements defined at the design stage, in the presence of the set of anticipated disturbances, without violating any process constraints. One of the basic requirements in performing operability analysis is that a model of the process relating the inputs and the disturbances to the outputs is available. Many of the process models can be described by the state-space representation

x3 ) f(x, u, d)

(1)

y ) g(x, u, d)

(2)

where x ∈ Rnx is the state vector, u ∈ Rnu is the input/ control vector, d ∈ Rnd is the disturbance vector, and y ∈ Rny is the output vector of the process. Here, x3 represents the time derivative of the state vector. It is implied in the above equations that all of the variables are functions of time. The two nonlinear maps f and g are of the dimensionalities Rnx+nu+nd f Rnx and Rnx+nu+nd f Rny, respectively. Additionally, we might have some constraints of the form

h1(x3 , x, y, u3 , u, d) ) 0

(3)

h2(x3 , x, y, u3 , u, d) e 0

(4)

representing the process, product, and safety specifications, including the bounds on the magnitudes and the rates of change of the inputs. These constraints might be applicable over the complete time history of the process (path constraints) or only at certain times (point constraints). The steady-state model of the process can be obtained by setting the time derivatives in eqs 1, 3, and 4 to zero. We discuss the steady-state and dynamic operability analysis techniques in the following two subsections. 2.1. Steady-State Operability. Before presenting the proposed modifications, we first review the concepts introduced by Vinson and Georgakis.2,3 The inputs of the process are each defined as being able to change over a range, the union of which is called the available input space (AIS). Subsequently, we can solve the model for the entire AIS to obtain the outputs that can be achieved; the collection of all of these output points is referred to as the achievable output space (AOS). Alternatively, suppose that we also specify a desired output space (DOS), which is the desired operating window for the process outputs. The set of input values required to reach the entire DOS can then be calculated from the inverse of the model given by eqs 1 and 2. The

4240

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

collection of all such input values is denoted as the desired input space (DIS). It is easy to see that the AOS is a function of u and d and the DIS is a function of y and d; this fact will be emphasized with subscripts and arguments. For example, the previously described AOS could be written as AOSu(dN) to imply that it is the AOS calculated by considering all of the points inside the AIS, denoted by the subscript u, when the disturbances are at their nominal values, dN. In the same way, the DIS described in the previous paragraph would be DISy(dN). The subscript y refers to all y ∈ DOS, and the argument shows the value of the disturbance used in the calculation. We do not use the subscripts and arguments if the context makes them evident. Having outlined the required spaces, we are ready to define the servo operability index in the output space as

s - OIy )

µ[AOSu(dN) ∩ DOS] µ[DOS]

(5)

Here, µ is a measure function that calculates the size of the corresponding space. For example, in two dimensions, it represents the area, and in three dimensions, it represents the volume. This index is more useful in analyzing the operability of existing plant designs, as it indicates how much of the desired process output region can be achieved with the available inputs. An index that is less than unity implies that our expectations of the process are higher than it can possibly deliver. For new plant designs or for the consideration of modifications to an existing plant, we can define a servo operability index in the input space as

s - OIu )

µ[AIS ∩ DISy(dN)] µ[DISy(dN)]

(6)

This index quantifies how much of the servo DIS is covered by the AIS. A value that is less than 1 indicates the need to increase the available ranges of some of the inputs. For linear systems, both of these operability index definitions give the same values, but there are some differences in the case of nonlinear systems. However, if one of the above indices were equal to 1, then the other would also be 1. To investigate the regulatory operability of the process, we need to specify the anticipated ranges of disturbances that define the expected disturbance space (EDS). In the steady-state case, the EDS might also reflect the uncertainties in some of the important model parameters employed in the design, such as kinetic constants, heats of reaction, heat-transfer coefficients, etc. The regulatory operability index is defined based on the inputs required to compensate for the effect of disturbances while maintaining the plant at its nominal set point, yN, as

r - OI )

µ[AIS ∩ DISd(yN)] µ[DISd(yN)]

theoretic union of DISd(y) for all y in DOS or the similar union of DISy(d) for all d in EDS. These two unions are identical to each other. Mathematically, we have

DIS )

∪ DISd(y) )

y∈DOS

∪ DISy(d)

d∈EDS

(8)

Then the overall operability of the process is defined as

OI )

µ[AIS ∩ DIS] µ[DIS]

(9)

Notice that all of these indices have best and worst values of 1 and 0, respectively. As they are ratios of similar quantities, these measures are dimensionless and do not suffer from scaling problems. Some variants of these indices might also be of interest. We could consider weighting different parts of the DOS and/or EDS depending on their relative importance. Another variant is that, in place of the DIS, we could use the nu-dimensional region constructed by calculating only the upper and lower bounds of each input variable in the DIS. Later in this section, we present an optimization-based method for calculating the bounds of the DIS. 2.1.1. Nonsquare Design Problem. As pointed out in the Introduction, Vinson and Georgakis2,3 and Subramanian and Georgakis9,10 considered systems with equal numbers of inputs and outputs. In a typical plant, however, one might have to decide on the appropriate values of more input variables than set-point-controlled outputs. Many of the output specifications would be inequalities, such as limits on the allowable pressure, temperature, and level in different process units. For instance, in the CSTR problem that will be studied later, the operator would require the plant to deliver a specific conversion (set-point controlled) but might not have such a strict restriction on the temperature or level in the reactor (set-interval controlled). It would be sufficient to keep the temperature and reactor level within some allowable limits. In this article, we propose an approach to deal with problems of this kind. The DIS calculation should be modified to account for the extra degrees of freedom. The strategy proposed here is to solve the steady-state design problem, that is, find the inputs, u*, for given ysp and d, as a constrained optimization problem minimizing a cost function subject to process and performance constraints. Mathematically stated

J(x*, u*) ) min J(x, u)

P1:

u

s.t.

(10)

f(x, u, d) ) 0 g1(x, u, d) ) ysp h1(x, y, u, d) ) 0 h2(x, y, u, d) e 0 ysp, d given

(7)

Our overall objective is to reject the expected disturbances and, at the same time, be able to reach all of the points in the DOS. With the understanding that DISd is a function of y and DISy is a function of d, the overall desired input space can be defined as the set

The objective function, J, would typically involve the cost associated with input use, soft penalties to discourage certain inputs, etc. The constraints of type g1 represent the equality output constraints, with ysp being the set point for those outputs. The constraints representing the set-interval-controlled outputs are included in h2. Note that u* in P1 is a function of ysp and d. To

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4241

calculate the overall DIS, this design problem has to be repeatedly solved, in principle for all ysp ∈ DOS and d ∈ EDS (see eq 8). Here, the DOS is specified only in terms of the set-point-controlled output variables. One should note that the servo DIS will have the same dimension as the DOS, which is smaller than the dimension of the AIS. The intersection calculations for obtaining the servo operability index defined in eq 6 can still be performed in the lower dimensions of the input space. 2.1.2. Bounding Box of DIS. The main strength of the approach developed by Vinson and Georgakis2,3 is its ability to represent the operational requirements of a plant in its entirety. For given servo and regulatory specifications, the overall DIS is explicitly calculated and presented. This information is of great value to process designers in evaluating a particular plant design, but as the dimensionality of the problem increases, the amount of data this approach generates becomes overwhelming. First, the method for performing the involved calculations is not always straightforward. Even the calculations were performed, it would be difficult to analyze and interpret such rich data. Therefore, a shortcut method to aid us in this analysis would be of great value. In the absence of the complete DIS, design decisions can be made on the basis of the extreme values of the input variables in the DIS. In other words, it might be of interest to identify the bounding box that encloses the DIS. Assuming that the inputs are independent of each other, we formulate the following optimization problem for finding the maximum demand in input uj

P2:

u j j ) max u/j (ysp, d) ysp∈DOS d∈EDS

(11)

Here, u/j (ysp, d) represents the jth element of u*(ysp, d), which is obtained by solving P1, and thus, P1 is nested in the above optimization formulation. This calculation has to be repeated for all of the inputs to obtain their individual upper limits. If the maximization operator in P2 is replaced by the minimization operator, we obtain the lower bounds on each input, uj. In addition to the input limits, the solutions of these optimization problems also identify the limiting combinations of outputs and disturbances that require these extreme values of the inputs. If the bounding box is found to be too large, the sensitivity of these input ranges to ysp and d can be used in an iterative effort to modify the DOS and the EDS to be realistic. 2.2. Dynamic Operability. The dynamic operability analysis technique utilized in this paper is the extension of the steady-state analysis presented in the previous subsection to the dynamic case. We briefly review this technique below. To perform the dynamic operability analysis effectively, we need quantitative measures of dynamic performance. Rise time, settling time, overshoot, and integral square error are well-known examples of such measures. Here, we use a measure similar to the settling time to quantify the dynamic performance. The dynamic operability measure is defined as follows:

Dynamic Operability Measure. The shortest time it would take a system to settle to the desired set point after a set-point change and/or a disturbance occurrence.

The operability measure is based on the idea that the time spent in driving the system back to the nominal operating point when a disturbance enters or in tracking a set-point change is directly linked to the cost of offspecification products. An estimate of this operability measure can be obtained using different types of controllers. However, a performance measure independent of the controller and capable of assessing the inherent limitations of the process is desirable. Time-optimal control suits these requirements very well. 2.2.1. Minimum-Time Optimal Control Problem. The minimum-time optimal control problem for continuoustime systems can be formulated as

∫0t dt

t/f (ysp, d) ) min u

f

(12)

s.t. x3 ) f(x, u, d); x(0) ) x0 y ) g(x, u, d) h1(y, x3 , x, u3 , u, d) ) 0 h2(y, x3 , x, u3 , u, d) e 0 x0, ysp, d given where t/f (ysp, d) is the minimum time necessary to respond to a change in the set point, ysp, and to a disturbance, d. The equations h1 and h2 represent the process constraints, including the bounds on the magnitudes and the rates of change of the input variables as well as the final-time constraints, which are typically set to ensure that the system reaches or returns to the set point, ysp, at the final time, tf, and stays there afterward, i.e., all state derivatives are equal to zero for all t g t/f . Note that we can set these constraints so that some of the outputs are controlled at the set point whereas others are controlled within an allowable range, as in the steady-state analysis. The bounds and rate-of-change constraints on the inputs define the equivalent AIS for the dynamic case, called the dynamic available input space (dAIS). For the optimal control problem in eq 12 to have a solution, it is necessary for the set-point changes to be tractable and the disturbances to be rejectable at steady state. In other words, the constraints must have at least one feasible solution at steady state; otherwise the optimal control problem has no solution. This is also equivalent to having a system that is steady-state operable, i.e., OI ) 1.0. The solution of the minimumtime control problem can be computationally intensive. Several numerical techniques exist in the literature for solving this optimal control problem. In most cases, it is transformed into a simpler problem, which is then solved iteratively for the minimum time and the corresponding control trajectories (e.g., see Gutman17 for a review of linear techniques). Here, we adopt a technique that converts the optimal control problem into a single nonlinear programming (NLP) problem. The details of this technique are given in Appendix B. 2.2.2. Dynamic Operating Spaces. Dynamic operating spaces are the extension of the steady-state operating spaces to the dynamic case. They include process constraints and specifications that need to be considered during transient operation. The first operating space, the dynamic available input space (dAIS) mentioned before, is defined as the set of values that the input variables and their combination can take. A typical dAIS

4242

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

S1 ) {(ysp, d)| ∀ysp ∈ DOS, ∀d ∈ EDS}

includes the constraints on the magnitudes and the rates of change of the inputs. The next operating space is the dynamic desired operating space (dDOpS), which is defined as the space formed by the combination of the DOS, EDS, and the desired response times. The mathematical definition of dDOpS is given by

The next operating space, referred to as S2, represents the ranges of set points and disturbances that can be achieved within tdf (ysp, d) or less. S2 is obtained by projecting the intersection of dDOpS and dAOpS onto the S1 and is defined by

dDOpS ) {(tf, ysp, d)| tf e tdf (ysp, d), ∀ysp ∈ DOS, ∀d ∈ EDS} (13)

S2 ) {(ysp, d)| t*f (ysp, d) e tdf (ysp, d), ∀ysp ∈ DOS, ∀d ∈ EDS} (16)

where tdf (ysp, d) represents the desired dynamic performance, or the maximum allowable response time, in tracking a set-point change, ysp, in the DOS and/or recovering from a disturbance, d, in the EDS. dDOpS can be easily formulated for the servo and regulatory problems separately by setting the variables that are not considered (disturbance and set point, respectively) to their nominal values. Note that, even though a lower bound for the response time tf is not specified, the response time has an implicit lower bound of zero. The dynamic achievable operating space (dAOpS) is defined as the operating space representing the dynamic performance that can be achieved by the system for a given choice of the dAIS, DOS, and EDS. Mathematically, dAOpS can be defined as

The dOI can be mathematically defined on the basis of these operating spaces as

dAOpS ) {(tf, ysp, d)| tf g t/f (ysp, d), ∀ysp ∈ DOS, ∀d ∈ EDS, u ∈ dAIS} (14) Because the response of a stable system to any (ysp, d) can take infinite time to reach the desired steady-state when u is changed to the corresponding steady-state value, dAOpS is unbounded at one end. However, a lower bound for the response time (i.e., an upper bound for the dynamic performance) can be obtained by using the measure of dynamic operability presented in the previous subsection. The calculation of this lower bound involves the repetitive solution of the minimum-time optimal control problem in eq 12 for each (ysp, d). Note that, if the system is not steady-state operable for some values of (ysp, d), then a solution to the minimum-time control problem, or t/f (ysp, d), does not exist at these values. In the following subsection, a dynamic operability index is presented for the operability assessment problem. The index exploits the concept of operating spaces to quantify the dynamic performance over the entire operating ranges. 2.2.3. Dynamic Operability Index. A comparison of the desired and achievable dynamic operating spaces (dDOpS and dAOpS) is used to define an index for dynamic operability. This index is referred as the dynamic operability index (dOI) and is defined as follows:

Dynamic Operability Index. The fraction of the operating ranges that can be achieved within the desired response time, t/f (ysp, d), given the available input ranges in dAIS. To mathematically define the operability index, two additional spaces are introduced. The first operating space, referred to as S1, is the space created by the combination of the set points in DOS and disturbances in the EDS

dOI )

µ(S2) µ(S1)

(15)

(17)

where µ represents a function calculating the size of the corresponding space. The dOI can take values varying between 0 and 1, representing the worst and the best performances, respectively. Note that the dOI can be calculated as a servo or regulatory dynamic operability index by selecting the appropriate dDOpS. Because the achievable response times are computed using an idealized controller, the dOI represents an upper bound for the achievable control performance of the process. That is, if the minimum-time controller fails to satisfy the performance requirements in the dDOpS (i.e., dOI < 1), then no feedback controller can satisfy them in practice. However, even if the performance requirements are met by the optimal controller (i.e., dOI ) 1), it does not necessarily guarantee that a realizable controller that will give the same performance exists. Further details about the dynamic operating spaces and dOI can be found in Uztu¨rk and Georgakis.11 3. Process Description In this paper, we study the interaction between design and control of single-CSTR and two-CSTRs-in-series systems with a first-order reaction of type A f B. The two configurations have the same feed flow rates and conversion specifications. Schematics of these systems are shown in Figure 1. For a liquid-phase exothermic reaction taking place in a jacketed CSTR, the mass and energy balances can be written, with assumptions of constant physical properties and complete mixing, as

Fcp

dV ) F0 - F dt

(18)

d(VCA) ) F0CA0 - FCA - VkCA dt

(19)

d(VT) ) FcpF0T0 - FcpFT + dt (-∆ H)VkCA - UA(T - Tc) (20)

FccpcVc

dTc ) FccpcFc(Tc0 - Tc) + UA(T - Tc) dt

(21)

The notation is defined in the Nomenclature section. The model equations are made dimensionless with appropriate transformations and are presented in Appendix A. In this investigation, we work with a system studied by Luyben.12 The parameter values are also listed in Appendix A.

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4243

Figure 1. Schematic of single-CSTR and two-CSTRs-in-series systems.

3.1. Variable Selection and Specification. In the steady-state analysis, the coolant flow rate, product flow rate, and reactor holdup are considered as the input variables. Even though the vessel volume of the reactor is a design variable, the volume of the reacting mixture is considered as an input variable in the steady-state operability calculations. This also helps us to possibly resize the vessel for given DOS and EDS specifications. The concentration of component A (or its conversion) in the reactor is selected as the set-point-controlled output variable, while the temperature is allowed to change within some bounds. The reactors are nominally designed to operate at an exit concentration of CA ) 0.05 lbmol/ft3, which corresponds to 95% conversion, a nominal feed flow rate of 100 ft3/h, and a feed temperature of 70 °F. The DOS is specified to be CA in the range [0.02, 0.2] lbmol/ft3. We let the temperature vary in the range 100 e T e 200 °F. A lower acceptable limit on the volume is fixed at 30% of the reactor volume. This avoids small-volume operations during low conversion demands. These bounds on CA, T, and V are used for both steady-state and dynamic operability calculations. The operating spaces DOS, EDS, and AIS used in this study can be represented mathematically as follows:

DOS ) {CA|CA ∈ [0.02, 0.2] lbmol/ft3}

| | | |

dF e 40 gal/min2 (or 19 250 ft3/h2) dt

(26)

dFc e 40 gal/min2 (or 19 250 ft3/h2) dt

(27)

In the case of two-reactors-in-series systems, the above definitions should be modified accordingly. The DOS needs to be defined only in terms of the concentration of component A in the second reactor, as this concentration determines the overall conversion. Because the EDS is defined in feed-related variables, it remains the same. The remaining eqs 24-27 should be considered separately for each of the reactors. Our preliminary dynamic calculations showed that the disturbances in the feed flow rate are more challenging than those in the feed temperature. Consequently, in the dynamic operability calculations, we considered the feed flow rate as the only disturbance to the process. Moreover, the effects of set-point changes and disturbances were evaluated independently. These results seem to represent the operability characteristics of the process well. If a different set of disturbances is thought of as important in some other application, they can be examined in a similar manner.

(22) 4. Single CSTR

EDS ) {(T0, F0)|T0 ∈ [50, 90] °F; F0 ∈ [50, 150] ft3/h} (23) AIS ) {(Fc, F, V)|Fc ∈ [0, 4FRc ] ft3/h; F ∈ [50, 150] ft3/h; V ∈ [0.3VR, VR] × ODF ft3} (24) Additionally, we have the constraint

100 °F e T e 200 °F

demands might change from time to time. The acronym ODF in eq 24 stands for the overdesign factor in the reactor volume, which will be obtained from the DIS calculations. This factor has a nominal value of 1. We discuss the issues related to the ODF in section 4.2. This is one of the main design decisions with which the operability analysis is expected to help us. It should be pointed out that, when the dynamic operability of this system is investigated, the holdup of the reactor (V) can not be used directly as an input variable. This is because it is not an independent variable. In fact, in the dynamic case, the reacting volume is an output variable that is to be controlled within the range given in eq 24. As a result, the set of input variables contains only the coolant flow rate and the product flow rate. Rate-of-change constraints on valve openings are also imposed in the dynamic calculations. It is assumed that all of the valves can open and close at a maximum rate of 40 gal/min2, which can be mathematically expressed as

(25)

The superscript R refers to reference or nominal values of the corresponding variables. Here, it is assumed that the maximum flow rate of the coolant is 4 times its nominal value. This was found to be sufficient in the steady-state operability calculations to be discussed later. The definition of the EDS is based on the anticipated disturbances in (T0, F0). The feed flow rate, F0, is considered as a disturbance because production

In this section, we present the design procedure, alternative design choices, steady-state operability considerations, and dynamic operability results for the single-CSTR system. 4.1. Design Procedure and Alternatives. The design of CSTR systems is well-treated in many reaction engineering textbooks. Here, we refer to some basic equations for the purpose of discussion. The design specifications include the feed flow rate, F0; the feed concentration, CA0; and the required conversion, x. The reactor exit concentration, CA, is related to the feed concentration and conversion through

CA ) CA0(1 - x)

(28)

Then, given the operating temperature, we can find the reactor volume using the steady-state material balance of component A as

4244

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Table 1. Three Single-CSTR Designsa,b T V Q Fc ∆T Qmax/Q cost Re(λ1) Re(λ2) Re(λ3)

D1

D2

D3

140.00 3800.00 2587.50 81.90 6.77 3.10 427.48 -0.025 -0.025 -18

160.00 1689.61 2512.50 63.88 11.29 2.91 258.07 -0.065 -0.065 -19

180.00 790.29 2437.50 53.12 18.17 2.67 160.79 -0.19 -0.19 -20

a F ) 100 ft3/h, C ) 0.05 lbmol/ft3. b F is in gpm, Q is in 103 A c Btu/h, and cost is in $1000 (using correlations in Luyben12). Refer to Nomenclature for units.

V)

F0(CA0 - CA) kCA

(29)

where k is the reaction rate constant at the reactor temperature. The reactor is sized assuming a heightto-diameter ratio of 2. The jacket volume is calculated taking a jacket thickness of 4 in. Finally, the nominal coolant flow rate, Fc, and jacket temperature, Tc, are calculated using the steady-state energy balances. The main design choice available in this problem is the operating temperature of the reactor. It can be seen from eq 29 that the volume of the reactor decreases as the operating temperature is increased. Therefore, it should be favorable to operate at a higher temperature from an economic point of view. If the cost of the cooling is also included, an intermediate optimal temperature might exist. In most applications, there are restrictions on the maximum allowable temperature based on factors such as safety, side reactions, material of construction, etc. In this work, we compare three single-reactor designs with nominal operating temperatures of 140, 160, and 180 °F, all with a nominal effluent composition of CA ) 0.05 lbmol/ft3 and a feed flow rate of 100 ft3/h. The nominal design conditions are given in Table 1. The controllability index of Luyben, Qmax/Q, is also shown in the table. These values are calculated using a maximum coolant flow rate of 4 times its nominal value. The real parts of the open-loop, nonzero eigenvalues of the linearized model, shown in the table, indicate that all such designs are open-loop stable. Note that the zero eigenvalue due to the overall material balance is not listed in the table. As discussed previously, one can see from Table 1 that the cost of the reactor decreases as the operating temperature is increased. Next, we present the steady-state and dynamic operability analyses of these designs. 4.2. Steady-State Operability Analysis. As described in section 2.1.2, each steady-state design requires the solution of the optimization problem P1 in eq 10. Here, we set the objective function J to minimize the reactor holdup. DISy(d) is obtained for a given d by solving P1 for the entire range of ysp ∈ DOS as defined in eq 22. Because the DOS is one-dimensional, the calculated DIS will also be a one-dimensional segment in the two-dimensional input space. Five different DISs are calculated and shown in Figure 2 for the 180 °F reactor (D3) for the nominal and four extreme values of the disturbances. Intermediate values of the disturbances lead to intermediate DIS segments. The union of all such DISs obtained for all of the disturbances in the EDS forms a two-dimensional space as shown in

Figure 2. Servo DIS corresponding to extreme values of disturbances in the EDS. The disturbance values (F0 (ft3/h), T0 (°F)) corresponding to different curves are nominal ) (100, 70), a ) (50, 50), b ) (150, 50), c ) (50, 90), and d ) (150, 90). The overall DIS is obtained as a union of many such DISs (see eq 8). The overall operability of the system is 68.2%. Table 2. Bounding-Box Vertices and Overdesign on Volume for Single-CSTR Designs V/VR a Fc/FRc a D1

D2

D3

0.30 0.30 0.30 0.39 0.30 0.30 0.30 0.90 0.63 0.30 0.30 1.90

0.31 1.67 0.59 0.84 0.32 1.88 1.06 1.00 0.38 2.31 0.60 1.17

CA

T

Tc

F0

T0

0.02 0.20 0.04 0.02 0.02 0.20 0.10 0.02 0.02 0.20 0.20 0.02

176.87 141.19 170.84 200.00 199.45 161.31 175.05 200.00 200.00 181.39 152.82 200.00

168.33 118.28 156.49 178.62 185.28 123.22 144.84 182.02 186.77 120.33 132.89 185.13

50.0 150.0 81.0 150.0 50.0 150.0 113.0 150.0 50.0 150.0 50.0 150.0

50.00 90.00 83.01 83.01 50.00 90.00 52.71 52.71 50.00 90.00 50.00 50.00

ODFb 5

8

128

a VR and FR are the reference values as given in Table 1. c ODF represents the required volume overdesign factor in terms of percentage of the nominal value.

b

Figure 2. It is worth mentioning here that, in the DIS calculations, whenever volumes greater than the nominal volume are demanded, the extra volume is thought of as being added through an increase in the height of the reactor. Next, the operability index is calculated by performing the required intersection calculation according to eq 9. In this case, we obtain an OI of 68.2%, implying that the design is steady-state inoperable. The DIS shown in Figure 2 reveals that, to guarantee steady-state operability, we would need 2.31 times the nominal coolant flow rate and 1.9 times the nominal volume. From the AIS, defined in eq 24 and also shown in the figure, one can infer that the available range of coolant flow rate is more than sufficient, but the available range of volume holdup is the limiting factor. Because it is desired for a process to have 100% steady-state operability, the reactor volume has to be overdesigned by 90%. These input bounds are also obtained using the bounding-box algorithm (P2) described in section 2.1.2. A direct search optimization technique (IMSL routine dbcpol) is employed, because of its better global convergence properties, for solving the outer optimization in P2. The results are shown in Table 2 for all three designs (D1, D2, and D3) considered here. Each row in the table represents a vertex of the bounding box.

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4245

Corresponding values of the state, output, and disturbance variables are also included in the table. The minimum and maximum of each of the input variables are emphasized with bold face. In the cases of designs D1 and D2, we observe that the maximum reactor volume is smaller than the nominal value. This is because the reactor temperature, T, is allowed to vary in a specified interval (see eq 25). The product flow rate, which is also an input variable, naturally takes on the value of the feed flow rate and thus is omitted from the table. In all cases, the available coolant flow range, between 0 and 4 times the nominal value, is found to be more than sufficient. The last column in the Table 2 records the overdesigns provided in the reactor volume. We note that dynamic operations demand input ranges larger than the steadystate requirements to achieve the desired dynamic performance in a finite time. Consequently, the reactor volume is taken to be 1.2 times the volume limit obtained in the bounding-box calculations. If the calculated volume of the reactor is smaller than 87.5% of the nominal volume, a marginal 5% overdesign is added to the nominal volume to accommodate disturbances that would temporarily increase the reactor holdup, such as an increase in the feed flow rate. For D1, the AIS is found to be more than sufficient to satisfy the steady-state operability demands for the given DOS and EDS. Still, a 5% overdesign is provided for reasons discussed above. Design D2 needs only an 8% overdesign (0.9 × 1.2 - 1), whereas reactor D3 required an overdesign of 128% (1.9 × 1.2 - 1). This might sound unreasonable, but in dimensional units, reactors D2 and D3 have the same total volume (∼1800 ft3). This is because the limiting case on the maximum volume corresponds to a reactor concentration of 0.02 lbmol/ft3 and a feed flow rate of 150 ft3/h, with both reactors operating at the maximum permissible temperature of 200 °F, as can be seen from the 8th and 12th rows of Table 2. This implies that, for given DOS and EDS, there is a minimum reactor size that is needed to have an operable process. As we briefly pointed out earlier, the extra volume needed for the overdesign is added by increasing the height of the reactor. Once the amount of overdesign is finalized, the reactor is resized so that its height-todiameter ratio is 2, as in the original designs. Here, we assume that this geometrical change only affects the coolant flow requirements because of changes in the heat-transfer area. This can be verified by performing a detailed operability study of each design as was reported by Russo and Bequette8,18 and Subramanian and Georgakis.9,10 Thus, the bounding-box calculations have to be repeated to find the new bounds of the inputs, in particular, of the coolant flow rate. These recalculated values are not reported here as they did not show significant differences. However, it should be noted that the normalizing coolant flow rate used in the dynamic calculations later corresponds to the modified design. 4.3. Dynamic Operability Analysis. Time-optimal control problems are formulated as nonlinear programming problems. In these calculations, the dimensionless nonlinear process model described in Appendix A is employed. The details of the solution procedure for the minimum-time calculations are presented in Appendix B. The minimum transition times that are needed for moving the reactor from the nominal concentration of

Figure 3. Minimum-transition-time plots for the three singlereactor designs. Transitions are from the nominal concentration of 0.05 lbmol/ft3.

Figure 4. Typical time-optimal output response of a single-CSTR design. This plot corresponds to design D3.

0.05 lbmol/ft3 to different concentrations in the DOS are calculated for the three designs under consideration. The results for the three reactors are compared in Figure 3. It can be seen that reactor D3 transits much faster than reactor D2, which transits faster than D1. Note that, in a small section of the high-conversion region (CA < 0.025 lbmol/ft3), reactor D3 becomes slower than designs D1 and D2. In this range, the reactor temperature in D3 reaches its maximum permissible limit of 200 °F, and thus, the required high conversion can be achieved only by increasing the reactor holdup. As volume build-up is a slower process, the transition takes longer. The observation that reactors operating at higher temperatures are faster than those operating at lower temperatures agrees with our intuition that reactors with smaller holdups have lower inertia (shorter residence times) and, thus, are able to move faster from one operating point to another. A typical time-optimal response for one of the transitions of D3 and the corresponding input moves are shown in Figures 4 and 5, respectively. In many instances, it was found that the final steady state reached by the reactor was unstable. Obviously, this means that a stabilizing controller should be in place to keep the system at its desired steady-state after the transition.

4246

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

These results clearly demonstrate that reactors with smaller holdups have inherently greater abilities for both servo and regulatory duties. This does not mean that one can design an arbitrarily small CSTR and expect it to perform better, even if the operating temperatures are acceptable. After a certain minimum volume, the available jacket heat-transfer area becomes a limiting factor. Moreover, it should be noted that minimum-time controller is an open-loop controller and any feedback controller will underperform. 5. Two CSTRs in Series

Figure 5. Typical time-optimal input trajectories corresponding to the output profiles shown in Figure 4.

The main motivation for selecting a system with two CSTRs in series over a single-CSTR system is the reduced capital cost. For a fixed reactor temperature, the captial cost of two reactors in series is much lower than the cost for a single large reactor. However, steadystate economics can conflict with the dynamic performance requirements. Therefore, the designer should not accept a process flow sheet solely on the basis of the steady-state economics. Both the steady-state and dynamic operability characteristics of processes should be examined before a final decision is made. In this section, we present the design procedure, design alternatives, and steady-state and dynamic operability analyses for systems with two reactors in series. 5.1. Design Procedure and Alternatives. As in the single-CSTR case, we start by specifying the desired feed flow rate, F0, and the desired exit conversion, x. Given the feed concentration CA0, we can calculate the concentration of component A in the effluent of the second reactor using the relationship

CA2 ) CA0(1 - x) Figure 6. Minimum-disturbance-rejection-time plots for the three single-reactor designs. Step disturbances in feed flow rate are considered.

The dynamic operability index defined in eq 17 can be calculated from these minimum-time plots. If the desired response time tdf is selected to be, say, 3 h for all ysp ∈ DOS (shown as a horizontal dashed line in Figure 3), then design D1 would have a dynamic operability index of 30%, which implies that only 30% of the set-point changes in the DOS can be achieved within the specified time of 3 h. Similarly, designs D2 and D3 would have dynamic operability indices of 60% and 99%, respectively. In this work, because we considered only one-dimensional dynamic operating spaces, time-optimal profiles can be graphically represented and compared. For this reason, we do not repeat the dOI calculations in the remainder of the article. Next, minimum-disturbance-rejection time calculations are performed for step disturbances in the feed flow rate for all three reactors. These results are presented in Figure 6. Interestingly, they show once again that reactor D3 is faster in rejecting the disturbances as well. This observation seems to be counterintuitive, as one would expect bulkier reactors to handle disturbances better than smaller ones. Actually, Zafiriou and Morari19 presented good arguments to show that set-point tracking and disturbance rejection need not be conflicting objectives, especially when step changes are considered. Similar discussions can also be found in Maciejowski.20

(30)

The steady-state material balances around each reactor lead to the following design equations:

F0CA0 F0 + V1k1

(31)

F0(CA1 - CA2) k2CA2

(32)

CA1 ) V2 )

Here, the temperature and volume in each of the reactors, as well as the concentration in the first reactor, are the unspecified variables. With these five unknowns and two equations (eqs 31 and 32), we are left with three degrees of freedom. We first choose to set the operating temperatures in the reactors, T1 and T2. One can possibly use the remaining degree of freedom to fix one of the following design variables V1, V2, or CA1. Instead, we select the volume ratio of the reactors as the last design specification so that we can observe its effect on the operability. This results in the additional design equation

V 2 ) R VV 1

(33)

where the volume ratio, RV, characterizes different design choices. For given T1, T2, and RV, eqs 31-33 can be rearranged to give the following quadratic equation in V1:

RVk1k2(1 - x)V12 + F0(1 - x)(k1 + RVk2)V1 F02x ) 0 (34)

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4247 Table 3. Steady-State Designs for Two CSTRs in Series with the Same Volumea,b T1 T2 V1 V2 CA1 Q1 Q2 Fc1 Fc2 ∆T1 ∆T2 Q1max/Q1 Q2max/Q2 cost Re(λ11) Re(λ12) Re(λ13) Re(λ21) Re(λ22) Re(λ23) a

VR10D1

VR10D2

VR10D3

VR10D4

VR10D5

VR10D6

VR10D7

VR10D8

VR10D9

140.00 140.00 694.43 694.43 0.22 2066.68 520.82 77.73 15.85 16.79 4.23 2.33 3.39 296.70 1.2 0.19 -21 -0.17 -0.17 -19

140.00 160.00 454.53 454.53 0.31 1820.81 691.69 72.34 16.77 19.63 7.46 2.17 3.20 227.88 2.1 0.054 -22 -0.27 -0.27 -19

140.00 180.00 294.72 294.72 0.40 1524.68 912.82 63.49 18.86 21.94 13.13 2.06 2.95 174.00 3 -0.1 -23 -0.5 -0.5 -20

160.00 140.00 454.53 454.53 0.16 2171.60 415.90 65.26 12.70 23.41 4.48 2.25 3.36 227.88 1.6 0.49 -22 -0.26 -0.26 -19

160.00 160.00 308.77 308.77 0.22 1991.68 520.82 64.06 12.60 27.78 7.26 2.08 3.22 179.12 3.1 0.19 -23 -0.4 -0.4 -20

160.00 180.00 207.75 207.75 0.30 1763.28 674.22 60.87 13.80 32.03 12.25 1.93 3.00 139.96 4.9 -0.049 -24 -0.68 -0.68 -20

180.00 140.00 294.72 294.72 0.12 2216.46 371.04 56.79 11.48 31.89 5.34 2.14 3.26 174.00 1.4 1.4 -22 -0.42 -0.42 -20

180.00 160.00 207.75 207.75 0.17 2087.06 425.44 57.94 10.35 37.91 7.73 1.97 3.18 139.96 4 0.59 -23 -0.58 -0.58 -21

180.00 180.00 144.42 144.42 0.22 1916.68 520.82 58.44 10.64 44.37 12.06 1.81 3.01 111.60 7 0.14 -25 -0.93 -0.93 -21

F ) 100 ft3/h, CA2 ) 0.05 lbmol/ft3, RV ) 1. b Fc’s are in gpm, Q’s are in 103 Btu/h, and cost is in $1000.

Table 4. Steady-State Designs for Two CSTRs in Series with a Volume Ratio of 2a,b T1 T2 V1 V2 CA1 Q1 Q2 Fc1 Fc2 ∆T1 ∆T2 Q1max/Q1 Q2max/Q2 cost Re(λ11) Re(λ12) Re(λ13) Re(λ21) Re(λ22) Re(λ23) a

VR20D1

VR20D2

VR20D3

VR20D4

VR20D5

VR20D6

VR20D7

VR20D8

VR20D9

140.00 140.00 484.43 968.86 0.29 1860.86 726.64 73.34 22.28 19.22 4.73 2.19 3.33 301.09 1.9 0.074 -22 -0.13 -0.13 -19

140.00 160.00 306.61 613.21 0.39 1553.15 959.35 64.44 23.55 21.77 8.47 2.07 3.12 226.47 2.9 -0.085 -23 -0.21 -0.21 -19

140.00 180.00 191.62 383.25 0.51 1205.42 1232.08 51.45 25.92 23.11 14.88 2.01 2.85 169.00 4 -0.31 -24 -0.41 -0.41 -20

160.00 140.00 327.30 654.60 0.21 2021.55 565.95 64.34 17.37 27.12 4.78 2.10 3.32 235.87 2.9 0.22 -23 -0.19 -0.19 -19

160.00 160.00 215.39 430.79 0.29 1785.86 726.64 61.27 17.76 31.67 8.12 1.95 3.15 181.77 4.7 -0.027 -24 -0.3 -0.3 -20

160.00 180.00 140.46 280.91 0.39 1499.46 938.04 54.92 19.54 35.36 13.94 1.84 2.90 139.28 6.9 -0.33 -25 -0.54 -0.54 -20

180.00 140.00 219.10 438.20 0.16 2108.85 478.65 57.79 14.80 36.97 5.29 1.99 3.26 183.71 3.6 0.68 -23 -0.3 -0.3 -20

180.00 160.00 149.30 298.60 0.22 1933.83 578.67 58.45 14.17 43.79 8.25 1.82 3.14 144.68 6.7 0.18 -25 -0.43 -0.43 -20

180.00 180.00 100.75 201.50 0.29 1710.86 726.64 57.40 15.06 50.35 13.47 1.69 2.93 113.25 11 -0.27 -27 -0.72 -0.72 -21

F ) 100 ft3/h, CA2 ) 0.05 lbmol/ft3, RV ) 2. b Fc’s are in gpm, Q’s are in 103 Btu/h, and cost is in $1000.

Once we solve for V1, then CA1 and V2 can be calculated from eqs 32 and 33, respectively. The nominal coolant flow rate and jacket temperature in each reactor can then be obtained from the steady-state energy balances. In this study, three nominal operating temperatures, namely, 140, 160, and 180 °F are considered for each of the reactors individually, resulting in nine different designs for a given RV. Additionally, we vary the volume ratio RV ∈ {1.0, 2.0, 0.5}. This helps us to investigate whether it is better to select two reactors with the same or different sizes. If different sizes are preferred, we wish to determine which reactor should be placed first in the series. The relevant details of these designs are given in Tables 3-5 organized based on the volume ratios. The real parts of the open-loop, nonzero eigenvalues of the linearized model are also included in these tables. Note that the process model has two additional zero eigenvalues, corresponding to the overall material balance for each of the two reactors. These are not listed in the tables. It can be seen from these tables that the first reactor in all of the designs is unstable at the nominal operating point, as evidenced by at least one positive eigenvalue. On closer inspection, we find that for some designs, for example, VR20D3 and VR20D6 in Table 4, where there is a good distribution of heat load between the two

reactors, only one of the eigenvalues is unstable, whereas in other cases, two of them are unstable. Taking this as an indication for systems that might have better controllability properties, we wanted to study the designs where the heat load is equally distributed between the two reactors. As most of the heat load is due to the heat of reaction, this equal heat distribution means the achievement of equal incremental conversions in each of the reactors. In Table 6, we present designs with equal heat loads for the same set of temperature combinations considered before. In this set, the volume ratio between the reactors varies over a broad range, depending on the operating temperatures of the reactors. One can observe that, in the cases where the second reactor is operating at 180 °F, the volume ratios are reasonable, and the cost of the designs are more desirable. Of the designs presented in Tables 3-6, eight were selected for a detailed operability study. Designs VR10D1, VR10D5, and VR10D9 from Table 3 are considered to act as our basis for comparison. In these cases, both reactors operate at the same temperatures of 140, 160, and 180 °F, respectively. The minimum-time results to be presented later show that, when both the reactors are operating at 180 °F, the system responds much more quickly. To understand the effect of volume ratio on this

4248

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Table 5. Steady-State Designs for Two CSTRs in Series with a Volume Ratio of 0.5a,b T1 T2 V1 V2 CA1 Q1 Q2 Fc1 Fc2 ∆T1 ∆T2 Q1max/Q1 Q2max/Q2 cost Re(λ11) Re(λ12) Re(λ13) Re(λ21) Re(λ22) Re(λ23) a

VR05D1

VR05D2

VR05D3

VR05D4

VR05D5

VR05D6

VR05D7

VR05D8

VR05D9

140.00 140.00 968.86 484.43 0.17 2224.18 363.32 80.16 10.98 14.48 3.75 2.47 3.45 301.09 0.5 0.5 -20 -0.23 -0.23 -19

140.00 160.00 654.60 327.30 0.23 2035.42 477.08 77.15 11.42 17.20 6.40 2.30 3.30 235.87 1.4 0.17 -21 -0.35 -0.35 -20

140.00 180.00 438.20 219.10 0.31 1797.36 640.14 71.73 12.97 19.85 11.22 2.16 3.06 183.71 2.1 0.042 -22 -0.63 -0.63 -20

160.00 140.00 613.21 306.61 0.13 2282.55 304.95 65.40 9.29 20.15 4.27 2.39 3.38 226.47 0.68 0.68 -21 -0.37 -0.37 -20

160.00 160.00 430.79 215.39 0.17 2149.18 363.32 65.18 8.70 24.01 6.44 2.22 3.29 181.77 1.8 0.43 -22 -0.53 -0.53 -21

160.00 180.00 298.60 149.30 0.23 1974.08 463.42 63.88 9.32 28.16 10.49 2.06 3.11 144.68 3.3 0.17 -23 -0.87 -0.87 -21

180.00 140.00 383.25 191.62 0.10 2293.78 293.72 55.78 9.13 27.70 5.63 2.28 3.22 169.00 0.8 0.8 -21 -0.61 -0.61 -21

180.00 160.00 280.91 140.46 0.13 2200.58 311.92 56.97 7.55 32.69 7.36 2.11 3.21 139.28 1.5 1.5 -22 -0.8 -0.8 -22

180.00 180.00 201.50 100.75 0.17 2074.18 363.32 58.02 7.32 38.46 10.69 1.95 3.10 113.25 4.2 0.54 -24 -1.2 -1.2 -22

F ) 100 ft3/h, CA2 ) 0.05 lbmol/ft3, RV ) 0.5. b Fc’s are in gpm, Q’s are in 103 Btu/h, and cost is in $1000.

Table 6. Steady-State Designs for Two CSTRs in Series with the Same Heat-Transfer Ratesa,b T1 T2 V1 V2 CA1 Q1 Q2 Fc1 Fc2 ∆T1 ∆T2 Q1max/Q1 Q2max/Q2 cost Re(λ11) Re(λ12) Re(λ13) Re(λ21) Re(λ22) Re(λ23) a

QR10D1

QR10D2

QR10D3

QR10D4

QR10D5

QR10D6

QR10D7

QR10D8

QR10D9

140.00 140.00 215.58 1725.00 0.48 1293.75 1293.75 55.00 40.29 22.93 5.73 2.02 3.21 333.04 3.7 -0.25 -24 -0.076 -0.076 -18

140.00 160.00 205.06 789.23 0.49 1256.25 1256.25 53.51 31.18 23.02 9.37 2.01 3.05 230.07 3.8 -0.27 -24 -0.17 -0.17 -19

140.00 180.00 195.06 379.55 0.51 1218.75 1218.75 51.99 25.62 23.09 14.81 2.01 2.85 169.13 3.9 -0.3 -24 -0.41 -0.41 -20

160.00 140.00 105.98 1625.00 0.46 1293.75 1293.75 48.68 40.43 36.81 5.96 1.80 3.19 297.90 8.4 -0.59 -26 -0.087 -0.087 -18

160.00 160.00 100.78 744.76 0.47 1256.25 1256.25 47.40 31.32 36.96 9.74 1.79 3.02 199.56 8.6 -0.65 -26 -0.19 -0.19 -19

160.00 180.00 95.86 358.75 0.48 1218.75 1218.75 46.09 25.78 37.08 15.38 1.79 2.82 141.56 8.9 -0.7 -26 -0.45 -0.45 -20

180.00 140.00 54.86 1525.00 0.43 1293.75 1293.75 48.95 40.60 57.10 6.22 1.56 3.16 272.66 17 -1.2 -29 -0.1 -0.1 -18

180.00 160.00 52.14 700.30 0.44 1256.25 1256.25 47.76 31.48 57.36 10.15 1.56 2.99 178.72 18 -1.3 -29 -0.21 -0.21 -19

180.00 180.00 49.57 337.95 0.46 1218.75 1218.75 46.50 25.95 57.55 16.01 1.56 2.78 123.42 19 -1.4 -29 -0.48 -0.48 -20

F ) 100 ft3/h, CA2 ) 0.05 lbmol/ft3, Q1 ) Q2. b Fc’s are in gpm, Q’s are in 103 Btu/h, and cost is in $1000.

system, designs VR20D9 from Table 4 and VR05D9 from Table 5 were also included. Finally, three designs, QR10D3, QR10D6, and QR10D9, from Table 6 were selected because the volume ratios and costs in these cases are more reasonable. 5.2. Steady-State Operability Analysis. Each steady-state design problem (P1) is solved with the objective function

J)

(

) (

)

V1 + V2 V1 V2 + -1 + -1 V1n + V2n V1n V2n

(35)

where V1 and V2 represent the reactor holdups of the first and the second reactors in a series, respectively. The additional subscript n indicates their nominal values, listed in Tables 3-6. The first term in this objective minimizes the total volume requirement of a design, while the last two terms penalize designs demanding excessive overdesign for either one of the reactors. These penalty terms were found to be necessary because a simple minimization of total volume always led to the same operable design, irrespective of the original design we started with, which is not desirable. As in the single-CSTR case, we first calculate the bounding box of the overall DIS so that each design can

Table 7. Overdesign in Volume for Selected Two-CSTRs-in-Series Designs VR10D1 VR10D5 VR10D9 VR20D9

ODF1a

ODF2

5 5 54 67

5 5 54 44

VR05D9 QR10D3 QR10D6 QR10D9

ODF1

ODF2

44 5 32 84

67 5 5 44

a ODF represents the required percentage volume overdesign i for reactor i.

be ensured of 100% steady-state operability. As described earlier, volume overdesigns are provided on top of the volumes obtained in the bounding-box calculations as a safety factor for operational flexibility during dynamic operations. These results are given in Table 7. Because we did not have the coolant flow rate in the objective, some lower-volume operations demanded excessively large coolant flow rates in the bounding-box search. This is remedied by constraining the coolant flow rate to be less than 3 times its nominal value. However, for the dynamic operability calculations, the AIS is assumed to be 4 times the nominal flow rate, as defined by eq 24. 5.3. Dynamic Operability Analysis. Minimumtime control calculations are performed for the set-point transitions from the nominal concentration to different points in the DOS for all eight selected designs. The

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4249

Figure 7. Minimum-transition times for three equal-volume twoCSTRs-in-series designs from Table 3.

Figure 8. Minimum-transition times for three two-CSTRs-inseries designs operating at 180 °F with different volume ratios.

results obtained for equal-volume (RV ) 1) designs are compared in Figure 7. This figure shows that the system with both reactors operating at 180 °F responds most rapidly. The results for the designs with different volume ratios, but with both reactors operating at 180 °F, are plotted in Figure 8. Interestingly, the figure shows that, for this system, the range of volume ratios examined does not have much effect on the response time. In other words, having the large reactor first or the small reactor first in the series does not have a significant effect on operability. However, certain caution should be exercised before this result is generalized further, as some additional calculations will indicate later. The minimum-transition-time plots for the designs with equal heat load distributions are shown in Figure 9. Here again, design QR10D9 with both the reactors operating at 180 °F responds most quickly. To make a comparison among the single-reactor and tworeactors-in-series systems, the minimum-transition-time curves of designs D3, VR10D9, and QR10D9 are plotted in Figure 10. This plot confirms that design QR10D9 has a significantly faster response when compared to the other two reactor designs. Moreover, it also reveals that both of the two-reactors-in-series systems are faster than a single reactor operating at the same temperature. Roughly speaking, QR10D9 is 5 times faster than VR10D9 and 10 times faster than D3. It should also be pointed out that the volume ratio for the two reactors

Figure 9. Minimum-transition times of three reactors-in-series systems, from Table 6, that are designed to share the heat duty equally.

Figure 10. Comparison of minimum-transition times of selected single-reactor and two-reactors-in-series systems. Two-CSTRs-inseries designs sharing equal heat duty show quicker responses.

in QR10D9 is equal to 6.8, indicating a more significant effect of reactor volume ratio than observed in Figure 8. Next, the minimum-disturbance-rejection time calculations for step changes in the feed flow rate for the eight selected two-CSTR systems are carried out. The results are presented in Figures 11-13. Among the equal volume designs depicted in Figure 11, design VR10D9 once again has the quickest response. This observation leads us to conclude that reactor systems operating at lower temperatures are slower than those operating at higher temperatures. Figure 12 shows that volume ratios in the range between 0.5 and 2 do not have a significant effect on the regulatory operability. Designs with equally distributed heat loads such as QR10D6 and QR10D9 are faster than most of the other designs, as in the servo case. This can be seen by comparing the response times given in Figure 13 with those given in Figures 11 and 12. This discussion points out that a good design criterion is even distribution of the heat load between the two reactors. The benefits could be substantial. 6. Conclusions This work proposed important methodological extensions to the operability framework of Vinson and

4250

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Figure 11. Minimum-disturbance-rejection times for three different equal-volume two-CSTRs-in-series designs with both reactors operating at the same temperature (from Table 3).

Figure 12. Minimum-disturbance-rejection times for three twoCSTRs-in-series designs operating at 180 °F with different volume ratios.

presented to identify the maximum and minimum requirements on each input, referred to as the bounding box of the desired input space, to ensure complete steady-state operability. These techniques were applied to single-CSTR and two-CSTRs-in-series systems to examine their steady-state and dynamic operability characteristics. Although the dynamic and control characteristics of CSTR systems have been a focus of many investigations in the past 30 years, only a few of these studies considered the effect of process design on control. Luyben12 and Jaisathaporn13 addressed this issue after selecting a particular control structure and controller tuning for the system. Although their approach is very useful, it does not elucidate the inherent dynamic characteristics of the system. In this work, following the recommendations of Uztu¨rk and Georgakis,11 we used a time-optimal control formulation in evaluating the true dynamic potential of the reactor designs. All of the designs studied were first made steady-state operable for the desired range of outputs, DOS, and the expected range of disturbances, EDS. Then, timeoptimal control calculations were performed to study both the servo and regulatory aspects of the inherent dynamic operability of the process. These studies showed that, in the single-reactor case, reactors operating at higher temperatures are able to transit more quickly to different set points and to reject disturbances more quickly than reactors operating at lower temperatures. In the two-CSTRs-in-series systems, the time-optimal control study showed that such systems are faster than single-reactor designs operating at the same temperature. When different designs of the two-reactor systems were compared, it was observed that the designs with the heat load evenly distributed between the reactors have only one unstable pole, instead of two as found in other cases. This was interpreted as indicating more favorable open-loop dynamics. The minimum-time calculations performed revealed that these designs are indeed faster and, thus, have significantly better dynamic operability characteristics. The results shown here are entirely based on openloop time-optimal control studies. The question of whether a feedback controller would deliver the predicted performance remains an open issue. In light of the results presented by Uztu¨rk and Georgakis11 for SISO systems, we might expect that an advanced modelbased controller, such as a model-predictive controller, would approach the performance of the time-optimal controller. The controller design problem for nonsquare systems of the kind studied in this article remains to be explored. Acknowledgment

Figure 13. Minimum-disturbance-rejection times for three reactors-in-series systems, from Table 6, that are designed to share the heat duty equally.

Georgakis2,3 and Uztu¨rk and Georgakis11 for dealing with nonsquare systems that arise when there are more input or control variables than set-point controlled outputs and when other outputs are set-interval controlled. This generalized framework helps us to investigate the operability of a wider class of design problems. Further, an optimization-based methodology is also

The financial support of the industrial sponsors of the Chemical Process Modelling and Control Research Center at Lehigh University and the Milestone Fellowship (to D.U.) is greatly appreciated. Nomenclature A ) heat-transfer area (ft2) As ) side heat-transfer area (ft2) Ab ) bottom heat-transfer area (ft2) cp ) heat capacity of A (Btu ft-3 °F-1) cpc ) heat capacity of coolant (Btu ft-3 °F-1) CA ) concentration of A (lbmol ft-3)

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4251 CA0 ) feed concentration of A (lbmol ft-3) d ) disturbance vector E ) activation energy (Btu lbmol-1) F ) flow rate (ft3 h-1) F0 ) feed flow rate (ft3 h-1) Fc ) coolant flow rate (ft3 h-1) J ) cost function k ) reaction rate constant (h-1) k0 ) Arrhenius constant (h-1) M ) molecular weight of A (lb lbmol-1) Q ) heat load (Btu) Qmax ) heat load at maximum coolant flow rate (Btu) RV ) ratio of volume in second reactor to that in first reactor RQ ) ratio of heat load in second reactor to that in first reactor tf ) response time T ) reactor temperature (°F) T0 ) feed temperature (°F) Tc ) jacket temperature (°F) Tc0 ) coolant feed temperature (°F) u ) input vector U ) overall heat-transfer coefficient (Btu h-1 °F-1 ft-2) V ) volume (ft3) Vc ) volume of the jacket (ft3) x ) state vector x ) conversion y ) output vector ysp ) set-point vector Greek Symbols F ) density of A (lb ft-3) Fc ) density of coolant (lb ft-3) λij ) jth eigenvalue of ith reactor (i is omitted for single reactors) (h-1) ∆H ) heat of reaction (Btu lbmol-1) ∆T ) temperature difference between jacket and reactor (°F) Superscripts/Subscripts i ) reactor number R ) reference value

A. CSTR Model Description. Mass and energy balances representing both single-reactor and reactorsin-series systems can be written as

dVi ) Fi-1 - Fi dt

(A1)

d(ViCAi) ) Fi-1CAi-1 - FiCAi - VikiCAi dt

(A2)

d(ViTi) ) FcpFi-1Ti-1 - FcpFiTi + dt (-∆H)VikiCAi - UAi(Ti - Tci) (A3)

FccpcVci

V R i Asi ) Asi VRi

(A5)

R represents the reference side heat-transfer where Asi area calculated at the reference volume, VRi , of the reactor holdup. Equations A1-A4 were made dimensionless with the following variable transformations:

x1i )

Vi

x2i )

VRi

qi ) φi ) β)

Fi FR

CAi

qci )

k0e-γVRi FR (-∆H)CRA

δbi )

CRA

FcpTR UAbi FcpF

R

Fci

R Fci

µi )

Fcp Fccpc

f(x3i) ) exp

τ)

FR

E RTR

ξ)

x4i )

TR

Ri )

R Fci

γ)

Ti - TR

x3i )

TR

tFR VR

Vi Vci

υi )

V1R VRi

R UAsi

δsi )

(

Tci - TR

FcpFR

)

γx3i 1 + x3i

The transformed equations are

dx1i ) υi(qi-1 - qi) dτ

(A6)

dx2i υi ) [qi-1(x2i-1 - x2i) - φi f (x3i)x1ix2i] (A7) dτ x1i dx3i υi ) [qi-1(x3i-1 - x3i) + βφi f(x3i)x1ix2i dτ x1i (δsix1i + δbi)(x3i - x4i)] (A8)

Appendices

Fcp

in the reactor (or proportional to the volume because a cylindrical geometry is assumed). This relationship is given by

dTci ) FccpcFci(Tc0 - Tci) + UAi(Ti - Tci) (A4) dt

where i represents the reactor number. The reaction rate constant, ki, is given by the Arrhenius equation ki ) k0e-E/RTi. The heat-transfer area, Ai, is calculated as the sum of the shell-side (Asi) and the bottom (Abi) heattransfer areas. Note that the side heat-transfer area changes with the reactor holdup. To account for this variation, we utilized the fact that the change in the side area is proportional to the height of the liquid level

dx4i ) υi[Riµiqci(x40 - x4i) + dτ ξµi(δsix1i + δbi)(x3i - x4i)] (A9) Here, x1i, x2i, x3i, and x4i are the state variables representing the normalized reactor holdup, concentration of reactant A, reactor temperature, and coolant temperature, respectively, in the ith reactor. The coolant flow rates and reactor volumes are normalized with the corresponding nominal values of each reactor. Therefore, in terms of the normalized variables, they have nominal values of unity. The system parameters are taken from Luyben12 and are listed in Table A1. B. Minimum-Time Optimal Control Problem. The minimum-time optimal control problem defined in eq 12 is solved using a NLP formulation. The first step in this procedure is to convert the infinite-dimensional problem into a finite-dimensional one. This is accomplished by discretizing the time interval t ∈ [0, tf] using N + 1 nodes such that the resulting N subintervals are of constant length T

T)

tf N

(A10)

4252

Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001

Table A1. CSTR Design Parameters parameter

value

units

E k140 °F CA0 U T0 Tc0 ∆H cp cpc M F Fc

30 000 0.5 1.0 300 70 70 -30 000 0.75 1.00 50.0 50.0 62.3

Btu lbmol-1 h-1 lbmol ft-3 Btu h-1 °F-1 ft-2 °F °F Btu lbmol-1 Btu lb-1 °F-1 Btu lb-1 °F-1 lb lbmol-1 lb ft-3 lb ft-3

Next, we assume that the control vector, u(t), has constant components on each of the subintervals. This assumption, in fact, means that a stepwise approximation to the actual control is being sought. Accordingly, the control vector is redefined as

u(t) ) u ˜ (k) t ∈ [kT, (k + 1)T], k ∈ [0, N - 1]

(A11)

and the unknown control moves are grouped to form ˜ T(N - 1)]T. the vector u ˜ (0: N - 1) ) [u ˜ T(0) ... u Further, it is assumed that it is sufficient to enforce the constraints h1(t) and h2(t) at discrete instants of time. Given an initial guess for u ˜ (0: N - 1), we can calculate the state, x(t), and the output, y(t), trajectories by integrating the system model. Then, the minimumtime optimal control problem can be cast as a NLP of the form

min

u ˜ (0: N-1)

tf

(A12)

s.t. h ˜ 1(0: N) ) 0

(A13)

h ˜ 2(0: N) e 0

(A14)

˜ T1 (t ) 0) ... h ˜ T1 (t ) NT)]T and h ˜ 2(0:N) where h ˜ 1(0:N) ) [h T T T ˜ 2 (t ) NT)] represent the constraints ) [h ˜ 2 (t ) 0) ... h h1(t) and h2(t) expressed at discrete instants of time. ˜ 2 are updated at each iteration via integration h ˜ 1 and h of the system model using the current u ˜ . Note that ˜ 2(t ) NT) represent the final-time h ˜ 1(t ) NT) and h constraints. We used the NPSOL21 optimization package to solve the NLP problems resulting from both the steady-state and dynamic operability analyses and the NAG routine d02eaf for the numerical integration of the nonlinear state-space model. All optimization problems were solved on an SGI Octane/SI R10000 workstation using Fortran 77. Each minimum-time calculation for a single CSTR takes about 20 min and for two reactors in series takes about 1 h on these computers. In the calculations reported in this study, we used N ) 25. For successful convergence, the need to provide a good initial guess cannot be overemphasized. Literature Cited (1) Morari, M.; Perkins, J. Design for operations. In Foundations of Computer-Aided Process Design; CACHE: New York, 1995.

(2) Vinson, D. R.; Georgakis, C. A new measure of process output controllability. In DYCOPS-5, 5th IFAC Symposium on Dynamics and Control of Process Systems; Georgakis, C., Ed.; Pergamon Press: Elmsford, NY, 1998. (3) Vinson, D. R.; Georgakis, C. A new measure of process output controllability. J. Process Control 2000, 10, 185-194. (4) Swaney, R. E.; Grossmann, I. E. An index for operational flexibility in chemical process design. Part I: Formulation and theory. AIChE J. 1985, 31, 621-630. (5) Dimitriadis, V.; Pistikopoulos, E. Flexibility analysis of dynamic systems. Ind. Eng. Chem. Res. 1995, 34, 4451-4462. (6) Bahri, P. A.; Bandoni, A.; Romagnoli, J. Effect of disturbances in optimizing control: Steady-state open-loop backoff problem. AIChE J. 1996, 42, 983-994. (7) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Integrated flexibility and controllability analysis in design of chemical processes. AIChE J. 1997, 43, 997-1015. (8) Russo, L. P.; Bequette, B. W. Impact of process design on the multiplicity behavior of a jacketted exothermic CSTR. AIChE J. 1995, 41, 135-147. (9) Subramanian, S.; Georgakis, C. Steady-state operabilty characteristics of reactors. Comput. Chem. Eng. 2000, 24, 15631568. (10) Subramanian, S.; Georgakis, C. Steady-state operabilty characteristics of idealized reactors. Chem. Eng. Sci. 2001, in press. (11) Uztu¨rk, D.; Georgakis, C. Inherent dynamic operability of processes: Definitions and analysis in the SISO case. Ind. Eng. Chem. Res. 2001, manuscript submitted. (12) Luyben, W. L. Trade-offs between design and control in chemical reactor systems. J. Process Control 1993, 3, 17-41. (13) Jaisathaporn, P. Trade-Offs Between Load Rejection and Switchability. M.S. Thesis, Lehigh University, Bethlehem, PA, 2000. (14) Douglas, J. M.; Denn, M. M. Optimal design and control by variational methods. Ind. Eng. Chem. 1965, 57, 18-31. (15) Douglas, J. M. The use of optimization theory to design simple multivariable control systems. Chem. Eng. Sci. 1965, 21, 519-532. (16) Vinson, D. R. A New Measure of Process Operability for the Improved Steady-State Design of Chemical Processes. Ph.D. Dissertation, Lehigh University, Bethlehem, PA, 2000. (17) Gutman, P. O. Controllers for Bilinear and Constrained Systems. Ph.D. Dissertation, Lund Institute of Technology, Lund, Sweden, 1982. (18) Russo, L. P.; Bequette, B. W. Operability of chemical reactors: Multiplicity behaviour of jacketted styrene polymerization reactor. Chem. Eng. Sci. 1998, 53, 27-45. (19) Zafiriou, E.; Morari, M. Setpoint tracking vs disturbance rejection for stable and unstable processes. In Proceedings of the American Control Conference; IEEE: Piscataway, NJ, 1987; pp 649-651. (20) Maciejowski, J. Multivariable Feedback Design; AddisonWesley: Reading, MA, 1989. (21) Gill, P. E.; Murray, W.; Saunders: M. A.; Wright, M. H. User’s Guide for NPSOL (Version 4.0): A FORTRAN Package for Nonlinear Programming; Systems Optimization Laboratory, Department of Operations Research, Stanford University: Stanford, CA, 1986.

Received for review December 13, 2001 Revised manuscript received March 19, 2001 Accepted March 21, 2001 IE001111+