Ind. Eng. Chem. Res. 1992,31, 2152-2164
2152
Williams, G.; Watts, D. C. Non-SymmetricalDielectric Relaxation Behavior Arising from a Simple Empirical Decay Function. Trans. Faraday SOC. 1970,66,80-85. Williams, M.L.; Landel, R. F.; Ferry, J. D. The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-Forming Liquids. J. Am. Chem. SOC.1955, 77, 3701-3707.
Zielinski, J. M.; Duda, J. L. Predicting Diffusion Coefficients in Polymer/Solvent Systems Using Free-VolumeTheory. AlChE J. 1992, 38 (3), 405-415.
Received for review February 3, 1992 Revised manuscript received May 22, 1992 Accepted June 16,1992
PROCESS ENGINEERING AND DESIGN Targeting Strategies for the Synthesis and Energy Integration of Nonisothermal Reactor Networks Subash Balakrishna and Lorenz T. Biegler* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
In this paper, a constructive approach for the synthesis of nonisothermal reador networks is presented based on a targeting methodology. Here, we determine a bound or target on a performance index for a given nonisothermal reacting system. The targeting model is based on mixing between different reacting environments and is formulated as a dynamic optimization problem, where the temperature, the feed distribution function, and an exit flow distribution function are the control profiles. In addition, the solution procedure is based on successive generation of reactor extensions that improve the target. Each reactor extension corresponds to the solution of a small nonlinear program; this procedure is repeated until no extensions that improve the objective function are generated. This technique ensures that only the simplest required model is solved. Having developed a general reactor targeting scheme, we also propose a new framework for integrating this scheme with an energy targeting approach. By discretizing the temperature profiles and allowing for heating or cooling profiles within the reactor, a combined representation for reactor and energy targeting is obtained. This model turns out to be a nondifferentiable nonlinear program, and novel smoothing techniques are presented for solution. Results on a typical process flowsheet optimization problem indicate significant increases in profit can be achieved by considering the reaction and energy synthesis schemes simultaneously. 1. Introduction
Most approaches to process synthesis have concentrated on the synthesis of chemical engineering subsystems. These include process subsystems such as heat exchanger networks, reactor networks, and separation systems. However, recently some attempts have been made toward integrating the synthesis of these subsystemsto generate systematic procedures for overall process synthesis. On one hand, many synthesis techniques based on hierarchical decomposition have been proposed (Douglas, 1985;Glavic et al., 1988). Though intuitive and simple to apply, they may not always lead to the optimal solution in the complete space of variables. On the other hand, rigorous mathematical programming techniques for overall synthesis are able to account for interactions between the different subsystems, but are usually limited by the size of the problems they can be applied to. In this paper, we first present a general mathematical programming technique for nonisothermal reactor network synthesis, and later develop an efficient scheme to integrate the concepts of reactor and energy network synthesis into a single framework. Even though the reaction system and the reactor design determine the character of the flowsheet, there does not seem to be a general approach for the synthesis of noni-
sothermal reactor networks. Most previous approaches to nonisothermal reactor network synthesis have been restricted to special forms of the reaction kinetics or objective functions. Aris (1960a,b,1961) applied dynamic programming to optimize the product yield or minimize reactor volumes in a variety of reactor trains. However, these were mainly limited to single reactions. Levenspiel (1962) suggests a geometric approach for single reactions, where the objective is usually the yield or residence time. Bhandarkar and Narasimhan (1969)proposed a method for the optimization of a train of adiabatic reactors with cold shot cooling, by successive optimization of single variable functions. Using variational calculus, Dyson and Horn (1967)propose a scheme for determining an optimal feed distribution to maximize yield for exothermic reversible reactions. The more recent publications include Chitra and Govind (19851,who classified different reaction mechanisms and presented a procedure for finding the optimal temperature profile in different reacting systems. However, their systems dealt mainly with yield maximization. In addition, Chitra and Govind (1985)provides a comprehensive review of the work done in nonisothermal reactor network synthesis. Recently, Glasser et al. (1990)extended their ideas on attainable region procedures (Glasser et al., 1987)to
0888-5885/92/2631-2152$03.00/00 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31, No. 9, 1992 2153 develop a novel geometric procedure for finding the entire region in concentration space for adiabatic reactions. On the basis of concepts of reaction and mixing, and by drawing alternate plug flow reactor (PFR) and continuous stirred tank reactor (CSTR) trajectories, they were able to map the entire attainable region in concentration space. However, since they follow a geometric approach, the techniques have dimensionality limitations. In this paper, we first develop a formulation for the synthesis of nonisothermal reactors, based on the constructive targeting approach proposed in Balakrishna and Biegler (1992). The model is a dynamic optimization problem, and simultaneous solution strategies are developed. The optimization algorithm resulta in the sequential solution of small nonlinear programming problems, where the solution to each program generates a component of the reador network. This provides a constructive technique for the target-based synthesis of nonisothermal reactor networks for any general objective function and process constraints. The synthesis technique is then illustrated on a few problems drawn from the literature. Reactors are invariably associated with heat effects that interact with the overall process. In the second part of this paper, we address the optimal heat integration of reactor networks based on minimizing the utility demands for the entire flowsheet. Today, energy network synthesis constitutes one of the most mature areas in proceas synthesis. Various techniques based on pinch technology (Linnhoff and Hindmarsh, 1983) have evolved over the years and have been applied successfully to a variety of energy integration problems. However, most of these considered fixed flow rates and temperatures of the heat exchanging streams. Recently, Duran and G r m a n n (1986) proposed a simple, yet rigorous framework for the simultaneous heat integration of streams with variable flow rates and temperatures in a flowsheet, based on utility bounding through pinch technology. Although pinch technology predicts strict bounds on utility consumption, overall sequential synthesis of the actual heat exchanger network using the pinch heuristic may lead to suboptimal solutions. Recent advances, especially the stagewise matching approach of Yee et al. (1990a,b), overcome this problem in a straightforward manner through a mixed-integer nonlinear programming formulation to model energy networks. However, this stagewise approach may be suboptimal for deriving pure utility bounds, because it assumes a fixed number of stages within which the heat exchange takes place. Despite these advances in techniques for energy integration, there has not been any work in the area of deriving reactor networks that are optimally energy integrated. Glavic et al. (1988) and Kravanja and Glavic (1989) performed a thermodynamicanalysis of reactors with fixed configuration and devised procedures to determine the relative placement of the reactor with respect to the pinch. They matched the process temperature-enthalpy profile with the reactor profile by controlling the reactor and process parameters to generate optimally heat integrated reactors. However, their technique was limited to either isothermal or adiabatic reactors, and they did not address the question of an optimal temperature trajectory within the reactor. With an efficient scheme for reactor network synthesis developed in the first part of the paper, we extend our synthesis approach to include the optimal energy integration of reactor networks. Here, we combine our reactor targeting formulation within a suitable energy targeting scheme,to derive a general framework for bounding energy costa while optimally synthesizing the reactor network
w
Figure 1. Segregated flow model.
within flowsheet constraints. Since our approach involves minimization of utility demands for the energy network, we apply the concepts of pinch technology within our reactor targeting framework. Solution strategies are presented and comparisons are drawn between sequential and simultaneous schemes for reactor and energy network synthesis. On a typical flowsheet example, significant differencesand savings in the reactor network are observed for simultaneous synthesis. 2. Model Formulation
The targeting model that we propose for nonisothermal reactors derives much of ita motivation from the constructive targeting methodology described in Balakrishna and Biegler (1992). In order to understand the generalization to nonisothermal reactors, we will briefly discuss the implications of that paper. We first start with the segregated flow (SF) model as shown in Figure 1, where a nonlinear programming problem to maximize the performance index in this reacting environment is described by the following formulation:
P1:
max J(Xerit,~) dXseg/dt = R(Xseg)
x m f ( t dt ) =1 Here, J is the performance index to be maximized within segregated flow. Xsegis the dimensionless concentration of molecules in the segregated environment, X o is the concentration of the feed that is entering the reactor denotes the reaction rate, t is the resinetwork, R(XSeg) dence time of a molecule, t , is the maximum residence time, f ( t )denotes the residence time distribution function, and 7 is the mean residence time. A segregated flow reactor is defined as one in which only molecules of the same age are well mixed within the reactor; all mixing between different molecules of varying ages can take place only at the reactor exit. A plug flow reactor (PFR) or a PFR with bypass falls into this class of reactors. For isothermal networks, such a reactor is described completely by the residence time distribution function, and this model can be reduced to a linear programming problem for yield- and selectivity-based objective functions. However, despite the restrictive mixing within segregated flow (SF),we observe that this model is sufficient for the synthesis of the reactor network in many cases. On the basis of these observations and some properties of attainable regions (Glasser et al., 1987; Hildebrandt et al., 1990), we can develop conditions for sufficiency of the region attainable from the segregated flow model. Moreover, even if these conditions are not satisfied, tests formulated as nonlinear programs can be used to evaluate the suitability and enhancement of the target obtained from the SF model. The solution to each
2154 Ind. Eng. Chem. Res., Vol. 31, No. 9, 1992
Reacting Segment I XI
To reador exif
Figure 2. Reactor extension from segregated flow model.
Figure 4. Reactor representation for discretized croea flow reactor model.
+
-I Figure 3. General cross flow reactor model.
of these nonlinear programs corresponds to reactor extensions from the region attainable to the SF model as shown in Figure 2. Successive convex hulls of these reactor extensions are formed in concentration space. Each point in this convex hull is feasible because creating the convex hull simply entails mixing the outlet concentrations from the reactors. Nonlinear programs (NLP’s) corresponding to reactor extensions are solved on each of these convex hulls, until we find no improvement in the objective function predicted by the solution of the NLP. These NLP’s usually tend to be very small because they are solved only on the variables corresponding to the reactor extension. The control variables in each iteration correspond to the residence time distribution function (which creates the convex hull on which the extensions are to be generated), and the parameters corresponding to the new reactor extension (Figure 2). In the synthesis of nonisothermal reactor networks, we adopt an approach in the same spirit of the ideas described above. However, one assumption in isothermal reactors is that mixing (or the creation of the convex hull) does not involve any significant costs. However, in the extension of these ideas to nonisothermal reactor network synthesis, there arises the cost of maintaining a temperature profile. One rather inexpensive technique for exothermic reactors is cold shot cooling. However, mixing may not always be optimal in the space of concentrations, even if it is desirable in terms of temperature manipulation. To address this, we consider a different reaction flow model as our basic targeting model, which can address temperature manipulation both by feed mixing as well as by external heating or cooling. The model consists of a cross flow reactor with a general exit flow distribution function. By construction, not only does this model allow the manipulation of reactor temperature by feed mixing, but it often precludes the need to check for PFR extensions, as will be described later. 3. Description of Targeting Model In this section, we develop the targeting model based on the ideas described above. Figure 3 shows a schematic of a cross flow reactor with side exists (CFR), which we choose as our basic targeting model. Here, Xois the concentration of the feed that is entering the reactor network. a is the independent variable denoting time as it progreaeea along the length of the reactor. T(a)denote8 the temperature as a function of the reactor length. We define f(a) da as the fraction of molecules in the reactor exit which left between points a and a + 6a of the reactor (an exit flow distribution function), and q(a) is the probability density function for a molecule entering the system at point a in the reactor. Thus the number of
molecules entering between points a and a da is given by q(a)Qo da, where Qo is the flow rate entering the reactor network. Implicit in the formulation is the assumption of instantaneous mixing between the feed and the mixture in the reactor. In addition, the formulation can easily be extended for variable-density systems, although we consider constant-density systems in this study. Clearly, at one extreme, when q(a) is zero throughout the reactor along with a general f(a), we have the equations for a segregated flow model. On the other extreme, when f(a) is a Dirac delta exactly at one point, and we have a general nonzero q(a),this model reduces to the Zwietering (1959) model of maximum mixedness. Glasser et al. (1990)have also considered a similar reactor without the side exits (the “differential”reactor) and concluded that, for single exothermic reversible reactions, this reactor may constitute part of the optimal network. Also, we defiie Q(a)as the flow of molecules within the reactor at point a. Based on this nomenclature, the mathematical model (P2) for maximizing the performance index in cross flow can be derived as shown below:
-dx - - R(T(a),X)+ *(X0 da Q(.)
imJ‘[q(a’)
- X(ff))
- f(a’)] da’ da = T
This formulation is an optimal control problem, where the control profiles are q(a),f(a),and T(a). The solution to this problem will give us a lower bound on the objective function for the nonisothermal reactor network along with the optimal temperature and mixing profiles. We now develop solution strategies for the above model based on orthogonal collocation (Cuthrelland Biegler, 1987)on finite elements. This type of discretization leads to a reactor network more practically achievable than the schematic shown in Figure 3 with a finite number of mixing points. It can be shown that if the finite elements (Aai) are chosen sufficientlysmall, the following formulation simply reduces to a numerical scheme for solving the system described above (seeAppendix). In the following section, we develop the discretized version of the model above, and Figure 4 shows the equivalent reactor representation of the discretized mathematical model.
Ind. Eng. Chem. Res., Vol. 31, No. 9,1992 2158
In the finite element discretization, the subscript i denotes the ith finite element, and the subscript j (or k) denotes the j t h (or kth) collocation point in any finite element. There are a total of N finite elements and K collocation points (i = 1,N,j = 1, K ) . The state variable X is approximated over each finite element by Lagrange interpolation basis functions &(a)) as K
X ( a ) = CXi&k(a) for aio I a I ai+l,o k-0
and
&(a)
=
ff a
i=0;k a i k
ail
ail
Substitution of this into P2 leads to the following nonlinear program, the solution of which gives us the optimal control variables at the collocation points. P3: max J(Xexitv7) +iJ",Ti
CXi&k'(aj) k
- R(Xij,Tij)Aai
0
X ( 0 ) = xo
(8) (9)
= CXikLk(ai+l,O)
(10)
4ixO + (1 - h ) X ( i - l ) e n d
(11)
x i end
xi,O =
= 1, K
k
CCfij 1 i J
OI4iIl (15) Here, Lk'(aj) are the derivatives of the Lagrange basis polynomials evaluated at aj, aj = (aij - aio)/Aai. Equation 8 repreaents the equations of orthogonal collocation applied to the differential equations at the collocation points. & correeponds to the ratio of the side inlet flow rate to the bulk flow rate within the reactor after mixing, and is an approximation to da)Qo/Q(a+)(4i = qiQo/Qi+l,l; Qij - ai0 is the length of each finite Cil{qG h i= ~ , and X , j correspond to the values of element r; f i j , c # J ~ Tij, the respective control and state variables at the jth collocation point in the ith finite element. The values for Xi, at an element are extrapolated to find the values x i a d at the end of that element in (10). Equation 11 represents the mass balance for the mixing between the feed and the reacting stream, while (12)and (13) are the integrals of P2 discretized using Gaussian qaudrature (see the Appendix). For convenience, the weights of Gauseian quadrature are absorbed into the respective discretized variables in the integration. Figure 4 is a schematic of the discretized reactor network. Thus the problem P2 can now be solved as a nonlinear program, to obtain the optimal set off, T,and I$ over each element. In this model, even though the temperature along the reactor is a control variable, part of the temperature manipulation can be readily accomplished by feed mixing if this is optimal for the reactor, and the remaining can be accomplished by the utilities. This will become more relevant when we consider energy-integrated reactor networks. Secondly, the cross flow reactor, by construction, mixea all available points on the reacting segments with the feed point and thus continuously checka for PFR extensions. So as long as the PFR trajectories are such that their nonconcavities can be enveloped from the feed point, one need not check for PFR extensions from the solution to this model.
The solution to this model provides a lower bound to the performance index of the reactor network. In the case of purely isothermal reactors, we can derive theoretical conditions of sufficiency for the segregated flow model (Balakrishna and Biegler, 19921,as well as check for favorable extensions from the segregated flow model, which are formulated as nonlinear programs. However, nonis o t h e d reactorsdo not lend themselves to easy analysis, because of the assumption of a completely arbitrary temperature profile. Nevertheless, by using some of the properties of attainable regions (Glasser et al., 1987),we can now develop techniques for extending the target provided by the CFR model. The constraints of P3 depict the feasible region for any achievable cross flow reactor. The convex combination of all the concentrations in thia region provides the entire region attainable by the cross flow reactor and mixing, which corresponds to the first candidate A, for the attainable region. Based on convex hull extensions illustrated in Figure 2,the following subproblem can be solved on A to check whether, starting from any point in A, a reactor can provide an extension that gives a point outside of A. Here, we consider a recycle reactor extension, since it includes the PFR and CSTR extensions as special cases. In the nonisothermal recycle reactor, we assume that the temperature is a control profile along the length of the plug flow section of the recycle reactor. The inlet temperature to the reactor which constitutes the extension will also follow a convex combination rule, if intermediate heating or cooling is not permitted. Recycle Reactor (RR) Extension: If J , > Jp3,then the recycle reactor provides an advantageous extension outside of A.
P4:
max Jrr(Xexit,7R) Xp3
= CCXijX,,(iJ
(16)
= R(X,,T,,)
(17)
i
dX,/dt
X,(t=O) =
l
+ XPI Re + 1
ReXexit
Here, J, is the value of the objective function at the exit of the recycle reactor. Xij is the convex combinator of all points available from the CFR model. Equation 16 describes the concentrations available from the cross flow reactor. The model equations for a recycle reactor starting from any point in 34 are described by (17)and (18).The variables X , and Re represent the concentrations in the recycle reactor extension and the recycle ratio respectively. X,, is the vector of exit concentrations for any reactor. Equation 19 gives the concentration at the exit of the recycle reactor. f , is a linear combinator of all the concentrations from the plug flow section of the recycle reactor (the indexes i and j correspond to collocation point j in element i for the discretized recycle reactor model), and 7 R is the residence time within the recycle reactor. Equations 20 and 21 are the typical normalization equations associated with the distribution functions for the recycle reactor extension. Here the vectors 1 and u are
2156 Ind. Eng. Chem. Res., Vol. 31,No. 9,1992
lower and upper bounds, respectively, on the exit concentration vector. Thus, if J, > Jp3,then the recycle reactor provides an advantageous extension over problem
P3. The next iteration consists of creating the new convex hull of concentrations, while including the concentrations obtained by this extension, and checking for favorable recycle reactor extensions from this new convex hull. At iteration P, this involves the solution of the following nonlinear programming problem: P5: max JP+’(u)
Q Sotve CFRModel
Is there a recycle
R.,&jJlowp),Tn(t)
dX,/dt
X,(t=O) =
= R(X,,TM)
(22)
+ &date Re+ 1
ReXexit
Find updated variables andobjecWe from
t 5 Form new convex hull of concentrations
P
P = P+l
I Here, Xmdel(p) is a constant vector and reflects the concentration at the reactor exit at iteration p in the models previously chosen. A convex combination of this with the cross flow region described by P3 gives the fresh feed point for the recycle reactor we are looking for, X,, then represents the concentration at the exit of the recycle reactor; and if J”+’) > J’p), the earlier model chosen is insufficient. The control profiles are f m d e l )I and T,, which are the linear combinators used to provize a convex candidate and the temperature profile in the recycle reactor, respectively. u is the set of d variables within the reactor targeting formulation. The solution to this NLP then gives the new reactor extension from the previous convex hull of concentrations. This procedure is repeated until no improvement in the objective function is observed between successive iterations. Figure 5 illustrates the flowchart for the synthesis procedure. It is easy to see that, with this constructive approach, the reactor network is synthesized readily. Also, due to the stagewise approach, only the simplest reactor model needed is solved. In other words, only the relevant part of the superstructure is generated, and thus many of the nonuniquenese problems with superstructure approaches are avoided. The example problems shown below illustrate the effectiveness of this approach for general nonisothermal systems. In fact, for many problems, the solution ie obtained through just one iteration of the flowchart in Figure 5, which is mainly because the variable temperature profile covers a large area for the attainable region. It can be shown that this algorithm will converge to the optimal solution whenever there are no extensions that generate a lower objective function and at the same time increase the region sufficiently so as to improve the objective beyond what we started with by solving a subsequent reactor extension. Results on practical problems shown below and in Balakriahna and Biegler (1992)suggest that this restriction does not seem to be limiting for the applicability of the proposed approach.
v,
4. Energy Integration of Reactor Networks Heat effects are associated with every reactor network, and energy integration involves matching heat loads between a set of hot streams and cold streams so as to
Figure 6. Flowchart for reactor synthesis.
minimize the cost of utilities for the network. In this section, we address the optimal integration of the heat effects within the reactor with the rest of the process. At the outset, we will first discuss the reactor flowsheet interactions. One can look at the reactor as the moat important unit of the chemical plant, because the downstream processing steps and the feed proceasing steps directly depend on the selective conversion of raw material to product. The main difference between stand-alone reactor network syntheais and reactor flowsheet integration is that the constrainta within the flowsheet allow a variable but usually bounded space of inlet concentrations for the reactor. Consider the case where the separation and feed processing unite are constrained so that the feed concentration to the reador is fixed. In this case,the constructive approach can be applied in a straightforward fashion. Since the inlet concentration to the reactor is known, the attainable region for the reactor in concentration space is independent of the flowsheet. Therefore, at each iteration of the synthesis algorithm only a problem of similar size needs to be solved. Each iteration correaponda to a reactor extension over the previous convex hull of concentrations, and since the previous hull is fixed, only the variables corresponding to the new extension are present in the optimization problem. These reactor extensions are of course generated while simultaneously satisfying the flowsheet constrainta. However, if the space of inlet concentrations is not predetermined, then at each iteration of the constructive approach, the problem size increases. This is because at each iteration, the previous reactors (or previous convex hulls) are not fued, as there is an infinite number of inlet concentrations dictated by the flowsheet constraints. Nevertheless, with either route to synthesis, the optimal flowsheet parameters can be determined by following the flowchart described in the previous section. In the following section, we discuse the development of a model for reactor flowsheet integration to ensure optimal energy management within the flowsheet. 4.1. Model Formulation. We consider two approaches for the integration of the reactor and energy network, the sequential and the simultaneous formulations. In the
Ind. Eng. Chem. Res., Vol. 31, No. 9, 1992 2167
4
Am Figure 7. Reacting segment for heat integration. Age
Figure 6. Optimal trajectory approximation.
conventional sequential approach, the reactor and separator schemes appear at a level higher than that of energy integration. In other words, once the "optimal" flowsheet parameters have been determined for the reactor target and the separation sptem, the reactor network is realized, and the heat exchanger network is derived in a straightforward manner. However, it is well-known that this approach is suboptimal with respect to the overall flowsheet. For the simultaneous approach, we consider reactor synthesis and energy integration at the same level. This approach is attractive because it considers the strong interaction between the chemical process and the heat exchanger network. However, this is not a trivial problem because the flow rates and the temperature of the process streams are not known in advance. Furthermore, the streams within the reactor cannot be classified as hot or cold streams a priori, because the optimal temperature trajectory within the reactor could be nonmonotonic. To address this, we discretize the temperature trajectory in the reactor and introduce the concept of candidate streams within the reactor network. Here, we approximate the optimal temperature trajectory within the reactor by a set of isothermal segments followed by the temperature change between these segments as shown in Figure 6. The solid curve represents the actual trajectory, while horizontal lines represent isothermal reacting segments and vertical lines represent the temperature changes in order to follow an optimal trajectory. The horizontal lines may correspond either to hot streams or cold streams depending on whether the reaction is exothermic or endothermic. The vertical sections may involve heating or cooling, and therefore we assume the presence of both heaters and coolers between the reacting menta Also, we call thew hot or cold streams candidate streams, because they may or may not be present in the optimal network. This will depend on the exit flow distribution in P3, which may require the inclusion of only some of the reacting segments and hence the corresponding temperature profiles. Figure 7 shows the reactor representation corresponding to the above approximation. This representation is similar to the discretized representation of the cross flow reactor shown in Figure 4, except for the additional heat exchangers. The subscript (superscript) i again represents the ith, finite element corresponding to the discretization. Tmi; corresponds to the temperature after mixing the reacting stream with the correspond the temperature enfeed. T h i b and Thiout tering and leaving the cooler, and Tihand Tid cornpond to the temperatures entering and leaving the heater. In an optimal network, it can be shown that only one of these two heat exchangers will be chosen. Aai corresponds to the length of the finite element, which may also be a variable in the optimization problem (subject to constraints on error control). Thus, the problem is well posed since we now know the hot and cold streams a priori, even if the
*To reactor exit
flow rates and the temperatures are not known. Also, while some amount of temperature control may be achieved by mixing, the remainder can be accomplished by the utilities or the heat flows within the network. Using the framework for reactor targeting from above, we now integrate this within a suitable energy-targeting framework. We assume here that utility costs will predominate and will be sufficient for the simultaneous total flowsheet synthesis for preliminary design. However, overall capital cost and area estimates (Kravanja and Grossmann, 1990)can also be included here if desired. As discussed earlier, several formalisms for energy targeting are available. Duran and Grossmann (1986) derived analytical expressions for minimum utility consumption as a function of flow rates and temperatures of the heat exchange streams. They showed that given a set of hot and cold streams, the minimum heating utility consumption is given by QH = max (zHp), where zHp is the difference between the heat sources and sinks above the pinch point for pinch candidate p . For hot and cold streams with inlet temperaturea given by Thinand t,'" and outlet temperatures Thoutand tcoUt, respectively, zHp(r) is given by wJmax(0, tFt- { P- ATJ - max(0, tcinzHpCy) = CEC
{ p- AT,)]] -
h€H
Wh[mU(O, Thin - p ]-
mU{o, - PI] (27) for p = 1,Np, where N is the total number of heat exchange streams. Here, & corresponds to all the candidate pinch points; these are given by the inlet temperatures for all hot streams and the inlet temperature added to AT, for the cold streams. w, and Wh are the heat capacity flows for the hot and cold streams (eqs 38 and 39). Also, y is the set of all variables in the reactor and energy network. The minimum cooling utility is given by a simple energy balance as Qc = QH + OCy), where My) is the difference in the heat content between the hot and the cold process streams, given by n(y) = xWh(Thin- Thout)- ~W,(TcoUt - Tp) h
C
On the basis of the above concepts for reactor and energy network synthesis, a unified target for simultaneous reactor-energy synthesis can now be formulated. We first classify the streams within the process into four categories. Let HR,CR be the set of hot and cold streams associated with the reactor network, and let Hp,Cpbe the set of hot and cold streams in the process flowsheet. Also, h E H = HR u Hp and c E C = CR Cp. If there are NE isothermal reacting segments, and if the reaction is exothermic, then there is a set of NE hot reacting streams from which the heat of reaction is to be removed in order to maintain a desired temperature in each segment. Also, between any two elements, there is a hot stream corresponding to the discretization shown in Figure 7. Hence, HR is a set of cardinality 2NE, while CR is a set of cardinality NE. For an endothermic reaction, CR and H R have cardinalities in the reverse order, since the reacting seg-
u
2168 Ind. Eng. Chem. Res., Vol. 31, No. 9,1992
menta now correspond to cold streams. Thus, there are always 3NE candidate streams, some of which may have zero heat content in the reactor network. Let Fh and F, denote the mass flow rates of the hot and cold streams, respectively. F(i)denotes the mass flow at the entry point of reacting segment i, and Fo is the total flow into the reactor. w h and w, are the heat capacity flows and Ah,,, Bh,,are the specific heat coefficients for temperature dependence. Here, the specific heats are assumed to be linear functions of the inlet temperature. w is the set of flowsheet parameters, and y constitutes the variables in the reactor and energy network, respectively, which include all variables in (281439) in the model below. Based on these assumptions, a unified target for simultaneous reactorenergy synthesis can be obtained by the solution of the following nonlinear programming problem: P6: *(w,QH,Qc) = ~ W Y -) CHQH - CCQC subject to
Cxi&k’(aj)- R(Xij,Tij)Aai= 0 k
j = 1, K (28)
X ( 0 ) = xo
xi,O =
(29)
9iXO + (1 - $ i ) X ( i - l ) e n d
(31)
= CCXijfij
(32)
Xerit
i l
Ci j f.. u = 7’
(33)
11
xfij =
v
F(0) =
F(ij) =
1
(34)
9Po
(35)
D i F ( i , o , - CfijFo ij
1
Ilei(t)ll -< 6 wh
=
(Ah
+ BhThin)Fh
Wh[Thin-
Thout]-
hEH
(36) (37)
w, = (A, + Bct>)Fc Q c = QH
0
(38) (39)
c W c [ t p t - t,’”] (40)
CEC
QH
2
ZH-YY)
(41)
h(wy)= 0
(42)
g(w,y) -< 0
(43)
where zHp, is given by the max expressions of (27). Here, TP corresponds to all the inlet temperatures for the hot streams and t,” + ATmfor all the cold streams. Thus (41) is of dimensionality 3NE + card(Hp) + card(Cp). The heats of reaction are directly accounted for by the heat capacity flow rates of the reacting streams as follows. If QR is the heat of reaction to be removed (or added, for endothermic reactions) to maintain isothermality in the reacting segment, the equivalent (FCJh or wh for this reacting stream is equated to QR, and we assume a 1 K temperature difference for this reaction stream. In the model above, (28)-(36) are discretization equations corresponding to the reactor network given by Figure 7. Note that in (33) we only evaluate an approximation f,instead of the actual residence time. This does not p a e a problem because the residence time only appears in the objective function, and it is expected that, for a general profit function, any minimizer for the actual residence time will also be a minimizer with respect to f . Equation 37
2
Figure 8. Nondifferentiable function: max(0,Z).
0
z
Figure 9. Hyperbolic approximation to 121.
ensures that the residuals (ei(t)for the differential equations) evaluated at a noncollocation point are within the tolerance 6. However, this equation is needed only if we allow variable lengths for the finite elements. Equations 38-41 along with (27) correspond to the heat integration equations, and (42) and (43) are the equations that bind the reactor energy integration variables (y) with the rest of the flowsheet variables (w). The solution to this problem will determine the optimal exit flow distribution and mixing profiles for the reactor network, while simultaneously minimizing the utility consumption for the entire flowsheet. 4.2. Hyperbolic Smoothing for the max Functions. It is clear that the above formulation corresponds to a nondifferentiable nonlinear program due to the presence of the max functions in (27). Here, we propose a simple smoothing technique based on a hyperbolic approximation to convert the above problem to a continuous nonlinear program. The max(0,Z)function has a nondifferentiability at the origin as shown in Figure 8. The main idea in our smoothing technique is to approximate the function 121, and then approximate max(0,Z) as f ( Z ) = max(0,Z) = lZ1/2 + Z/2
Also, the equation to a standard hyperbola is given as y2/P - x2/a2= 1, where a and b scale with the minor and major axis, y and x are the ordinate and abscissa, and fb/a correspond to the slopes of the asymptotes. For the function 121the slope of the asymptotes to the hyperbola will be f 1 and, if a is sufficiently small, the hyperbola will approximate 1213 shown in Figure 9. Thus, our approximation to the max function can be given as follows:
Besides having a simple form, this approximation also provides a single function representation over the entire domain (-m < 2 < m), unlike previous quadratic or exponential approximations. With values of t = 0.01, we were able to obtain a reasonably good approximation to the max function while simultaneously allowing easier solution of problem P6. 4.3. Extensions from This Targeting Model. The solutions of the nonlinear optimization problem P6 gives us a lower bound on the objective function for the flowsheet. However, this cross flow model may not be suffi-
Ind. Eng. Chem. Res., Vol. 31, No. 9,1992 2169 cient for the network, and we need to check if there are reactor extensions that improve our objective function beyond those available from the cross flow reactor. Here, we check for CSTR extensions from the convex hull of the cross flow reactor model, much in the same spirit as the illustration in Figure 2,except that all the flowsheet constraints are included in each iteration. In earlier work (Balakrishna and Biegler, 1992),we found that these extensions give quite good targets and yield much smaller problems. A CSTR extension to the convex hull of the crOss flow reactor constitutes the addition of the following instead of CP: constraints to P6 in order to maximize I
max CP(2)(wy(2),QH,Qc) = J ( U , ~ ( ~cHQH ) ) - - cCQc
(45) Xcstr = Xerit + R(X,,,T,,) T,, 7 2 0, X,tr 2 0 Here, X,, corresponds to the new reactor extension and y(z)is the set of new variables in the reactor and energy network, respectively. Beaides the set y in P6,this includes the variables corresponding to the new CSTR extension, namely, X,,, T,, T, and three more candidate streams for heat exchange. This is because of the two heat exchangers that will either cool or heat the feed to the CSTR (only one of these will exist in the optimal network), and one additional exchanger within the CSTR to maintain a 1 CP*, we have a reactor exdesired temperature. If tension that improves the objective function. The next step consists of creating the new convex hull of concentrations, and checking for extensions that improve our objective function within the flowsheet constraints. As in section 3,we continue this procedure until there are no extensions that improve the objective function. In the following section, we demonstrate the approaches presented above on a few example problems. All examples were solved using the GAMS modeling system (Brooke et al., 1988). 5. Example Problems
In thia section, four examples are presented to illustrate our targeting approach. The first set of three examples involve nonisothermal reaction mechanisms, drawn from the literature. The last example involves a detailed analpis of simultaneous reactor and energy network synthesis for a typical flowsheet. Example 1. The Van de Vusse reaction involves four species for which the objective is the maximization of the yield of intermediate species B, given a feed of pure A. The reaction network is given by kl
A - B - C
I
I
1.o
0.5
1.5
time(a )secs
Figure 10. Temperature profile in Van de Vusse reaction.
An upper bound to 500 K was used for the temperature and a bound of 2.5 s was placed on the residence time as in Achenie (1988). Following the algorithm described in Figure 5,we find that an optimal yield of B = 0.786 can be obtained in a recycle reactor (recycle ratio = 0.863,PFR section space time = 0.22 s) with a temperature profile in the plug flow section as shown in Figure 10. In comparison, Achenie (1988)obtained a maximum target yield of 0.1797 in a reactor network with a residence time of 2.0 8.
Example 2. This example is a case of parallel reactions taken from Kramers and Westerterp (1963). The reaction involves four species with the following reaction diagram:
The rate vector for components A, B, C, and D, respectively, is given by
where
ki= kio exp(-Ei/RT),
i = 1, ..., 3
and klo = 4.86 X 106 s-',
k2o = 2.2 X 10' k30
k2
k31
= 1.578 X 10"
5-l
El = 12 kcal/mol, E2= 7 kcal/mol, E3 = 18 kcal/mol
D
The reaction rate vector for components A, B, C, and D, respectively, is given by
R(C) = [-klCA - k3CA2,
0
klCA
- k2CB, k2CB,
k3CA2]
where
ki = kio exp(-Ei/RT),
i = 1, ..., 3
and klo= 5.4 X lo9 h-l, k20 = 1.6 X 10l2 h-l, k30 = 3.6 X lo5 mol-' h-', CAO= 1.0 mol/L
El= 15.84 kcal/mol, E2= 23.760 kcal/mol, E3 = 7.920 kcal/mol
Also, E2IEl IEa,and the objective is the masimization of species B, with a constraint on the residence time of 1.34 s. With these constraints, our constructive technique predicts a maximum yield XB = 0.3914. The optimal reactor network is a PFR with a residence time of 1.34 s with the temperature profile given in Figure 11. However, for unconstrained residence time, since the activation energy for the desired reaction is in between the activation energies for the undesirable reactions, we would expect that there is an optimal isothermal temperature to maximize selectivity. This was observed analytically by Kramers and Westerterp (1963). Example 3. Here, we maximize the conversion in the oxidation of sulfur dioxide, which has been investigated
2160 Ind. Eng. Chem. Res., Vol. 31, No. 9, 1992 630
Unnactcd A(!39%1
I
v
620
Q
RE4CTOR
580 570
Hl-Hl3 c2-c7
J
I 0
2
1
BCD ~
X
L
al
Figure 13. Flowsheet for reactor-energy network synthesis.
time (Q ) secs
Hence the cost of heating cannot be readily compared to their results, and thus we choose the extent of reaction as the objective function. For this reaction, we observed that for a constraint on residence time of 0.25 s, the maximum reaction extent of 2.42 is obtained in a PFR with the temperature profile shown in Figure 12A. However, if the constraint on the residence time is removed, the extent of reaction (as defined by g) asymptopically approaches the upper bound of 2.5 for a suitable temperature profile in a PFR with a sufficiently large residence time. Also, for a residence time bound of 2.2 s, the temperature profile as a function of time is shown in Figure 12B and delivers an extent of reaction of 2.48. Example 4. This example illustrates the simultaneous synthesis of reactor and energy networks. Here, we consider a reaction mechanism in the Van de Vusse form, though with different kinetic expressions. The integrated flowsheet correspondingto the synthesis problem is shown in Figure 13. The reactions involved in this flowsheet are as follows:
Figure 11. Temperature profile for example 2. 1100
A loo0 -
7
I
-I
I
d
0
I
200
100
300
time (a x1OOO) secs
1100
B
ma
1
I
4
k2
kl
600
-
A - 0 - C
"i
h a m .. .. .
500 0
1000 time (Q
2000
D
3000
where
x1000) secs
Figure 12. (A) Temperature profile for example 3. (B)Temperature profile for example 3 (large residence times).
klo = 8.86 X lo6 h-',
by Lee and Aris (1963). The reaction and the kinetics are described as follows:
El = 15.00 kcal/mol, Ez = 22.70 kcal/mol, E3 = 6.920 kcal/mol
1
SOz + 50' = SO3
k20 = 9.7
AHA+, = -0.4802 kcal/mol,
-0.918 kcal/mol,
R = 3.6 X 50 112.5 - g)0.5(3.46- 0.5gl 1 + 0.311t (32.01 - 0.5g)',5
-
1
Here, g is defined as the number of moles of SO3 formed J,where per unit m m of mixture, t is defined as (T- To)/ Tis the temperature, Tois 310 K (fresh feed temperature), and J = 96.5 K kg/mol. R, the rate of reaction is defined as the kilomoles of SO3 produced per hour per kilogram of catalyst. The extent of reaction or the moles of SO3 formed per unit mass of mixture is limited by the inlet mass flow of SO2,which is fixed at 2.5 mol of SO2per unit mass of mixture. Lee and Aria looked at the maximization of an objective fundion based on the value of the product stream, catalyst, and preheating costa. However, they assumed adiabatic reactor sections, with cold shot cooling. In our procedure, however, we do not enforce the adiabatic restriction.
X
k30
lo9 h-l, = 9.83 X lo3 mol-' h-'
AHB-c = = -0,792 kcal/mol
mA-D
The feed to the plant consists of pure A. This is mixed with the recycle gas stream consisting of almost pure A, and preheated (Cl) before entering the reactor. After reaction, the mixture of A, B, C, and D passes through an aftercooler prior to distillation. In the first column, A is recovered and recycled, while in the second column, the desired product B is separated from CD which is used as fuel. The distillation columns are assumed to operate with a constant temperature difference between reboiler and condenser temperatures (Andrecovich and Westerberg, 1985). The reflux ratios are fixed and the column temperatures are functions of the pressure in the column, which is variable so that efficient heat integration can be attained between the distillation columns and the rest of the process. The reactor is represented by the cross flow discretization shown in Figure 7. NE, the total number of reacting segments, is chosen to be 7, and the segment lengths Aa,are assumed constant. Since the reaction is exothermic, this corresponds to 14 hot streams and 7 cold streams. Thus the streams in the reactor may be enumerated as hot streams H1-H13 (2NE - 1 since the entry
Ind. Eng. Chem. Res., Vol. 31, No. 9, 1992 2161 SteFm
Hfl
4
SO0
C 480 460
0.0
0.2
0.4
0.6
0.0
0.2
0.4
0.6
Figure 14. Reactor temperature profiles for example 4. Table I
overall profit, $/yr overall conversion, % energy consumption, BTU/h hot utility load cold utility load
sequential simultaneous 38.98 X lo6 74.02 X lo6 49.6 61.55 3.101 X lo5 252.2 X lo6
2.801 X lo6 168.5 X lo6
point is fixed to be a preheater), and cold streams C 1 4 7 . The streams H15-H16, and C8, C9 correspond to the condensers and reboilers of the distillation column. As described in P6, the specific heats are assumed to be linear with the inlet temperatures. The objective function here is the total profit for the plant and is given by
J = 1 . 7 F ~+ 0.8FCD - 6.95 X 10-5~’F00.4566FA(1+ O . O l ( P H I 5 - 320)) - 0.7(FB+ F C D ) O.U;l,o- 0.007Qc - 0.08QH In this expression, F B and FcD represent the production rates of B and CD, respectively. FAois the amount of fresh feed. The third term corresponds to the reactor capital cost, while the fourth and the fifth terms correspond to the capital cost of the distillation columns. The operating costs of the columns are directly incorporated into the energy network in terms of condenser and reboiler heat loads. We assume that the cost of the reactor can be described by the total residence time and is independent of the type of reactor. This can be justified on the grounds that the capital cost of the reactor itself is usually orders of magnitude smaller then the operating costa and the capital costs of the downstream processing steps. A target production rate of 960 OOO lb/day is assumed for the desired product B. Here, we consider two alternatives. Firstly, we consider the sequential approach, where we optimize the reactor network to find an optimal temperature profile, and then integrate the maintenance of this optimal profile with the energy flows in the rest of the flowsheet. In the second case, we solve the above problem with the formulation proposed in P6. Here, the reactor and the energy network are optimized simultaneously. Table I gives a comparison between the results for sequential and simultaneous modes to synthesis. The optimal reactor in both cases, has a residence time of 0.59 s. However, the temperature profiles and the flow rates in the two cases are different. In the simultaneous case, of the 20 candidate streams, only 12 streams are active in the optimal network. This is because the strictly falling temperature profiles avoid the use of any cold streams ( C 2 4 7 ) within the reactor network. Also, no mixing was predicted in the solution, so cold shot cooling was not used at all. However, this is directly influenced by the ratio of the raw material to energy costs. For small ratios, even if mixing is not optimal in concentration space, the energy costa may drive the use of cold shots, to reduce utility consumption. It is clear that even though the temperature profiles in the two cases are not markedly different (Figure 14), the signifcant increases in profit are
Reacting Segment (Aa
)
-
2162 Ind. Eng. Chem. Res., Vol. 31, No. 9,1992 Table 11. Stream Data for Network (Temperature [K], Heat Capacity Flow Rate [lo' BTU/K/h] sequential simultaneous stream Tin To,, heat cap. flow rate stream Tin Tout H1 551.5 523.8 4.4 H1 540.1 512.3 4.2 H2 512.3 510.6 498.6 523.8 H2 487.7 4.1 H3 498.6 500.6 510.6 H3 4.0 H4 487.7 494.2 480.3 500.6 H4 490.0 475.2 3.9 H5 480.3 494.2 H5 6781.2 H6 540.1 550.5 539.1 551.5 H6 H7 512.3 522.8 511.3 2967.8 523.8 H7 HB 497.6 4150.6 498.6 509.6 510.6 H8 H9 4341.4 487.7 499.6 486.7 500.6 H9 H10 479.3 4048.8 480.3 494.2 493.2 H10 H11 415.2 489.4 474.2 2955.6 490.4 H11 3.9 H14 475.2 490.4 320.0 321.0 H14 H15 321.0 6.9 322.0 320.0 323.0 H15 854.8 H16 349.0 350.0 349.0 350 H16 c1 540.1 297.1 3.1 551.5 294.3 c1 C8 192.4 352.0 351.0 354.0 353.0 C8 c9 854.8 391.0 390.0 391.0 390.0 c9
either with oscillatory initial conditions or a smooth exponential exit distribution. Since the CFR model includes the Zwietering model, the CFR model can exhibit multiple steady states. In this paper, our algorithm proceeds through the discretization of the CFR model and does not seek to find all these solutions. This does not pose a major concern mainly for two reasons. Firstly, from the target to the CFR model, we derive recycle reactor extensions (which can be either PFR or CSTR extensions),where the multiple solutions that might be lost in the CFR model may be included in the attainable region after the extensions have been created. Secondly, since the problem is optimizationdriven (performance index maximiation), we expect to generate the most rewarding among these solutions, and use this as our candidate for the new reactor extension. Furthermore, since we allow a completely variable temperature profile (unlike typically adiabatic situations that have been analyzed for multiple steady states), the reaction curve itself will seek to proceed in a direction to maximize the objective. The emphasis, therefore, is mainly on generating the solution that maximizea the objective, rather thanenumerating the different solutions possible for the system. It must be noted, however, that most reactor synthesis problems can be highly nonlinear and that global optimality usually cannot be guaranteed. The main strength of this approach compared to previous approaches is that while postulating a suitably rich reapresentation for a reactor network, only the simplest model that is needed for the particular reaction system is solved. Also, with a target-based approach, nonuniqueness problems associated with superstructurebased approaches can be avoided, as only the simplest relevant part of the superstructure is generated. Examples 1-3 illustrate the successful application of these ideas on nonisothermal reaction kinetics. Particularly in example 3, we are able to show that the performance index actually reaches ita upper bound aymptotically as the residence time is sufficiently increased with a suitable temperature profile. In all these examples, the solution is obtained through just one iteration of the flowchartof Figure 5 with better results than those reported in the literature. In the second part of this paper, we have developed a novel scheme for deriving reactor networks while ensuring optimal energy management within the entire flowsheet. Here again, we describe two approaches based on sequential and simultaneous synthesis. If we look at the reactor as the most active unit of a chemical plant, then in simultaneous synthesis we have an infinite number of inlet concentrations for the reactor based on the con-
heat cap. flow rata 3.9 3.7 3.5 3.5 3.4 4653.2 2030.0 2831.1 2921.2 2691.3 1912.2 3.4 11.0 854.8 2.8 309.4 854.8
strainta in the flowsheet. The application of our targeting procedure is still essentiallythe same in spirit, except that the flowsheet constraints are embedded in each iteration and thus the complexity of the problem increases, as we described before. Through the discretizationof the reactor network, we are able to provide a classification of the streams within the reactor, while simultaneouslyreducing the system of equations to an algebraicone. The resulting model turns out to be a nondifferentiable nonlinear program, and a new smoothing technique is presented for easy solution. Example 4 illustrates the effectiveness of a scheme for the simultaneous synthesis of reactor and energy networks. In this case study, cold shot cooling was not predicted at the optimal solution. For any reaction system where mixing is nonoptimal, it could be conjectured that the use of cold shota for temperature control is mainly influenced by a high ratio of energy to raw material costa. Our next goal is to extend these ideas to include reaction-separation systems. Here, we are investigating techniques to analyze costa associated with separation on a thermodynamic basis. Preliminary results show that, in some casea, the conventional approach of a reaction system followed by separation can be suboptimal. These issues will be addressed in a future publication.
Acknowledgment This work was supported by the Engineering Design Research Center, an NSF supported center at Carnegie Mellon University.
Nomenclature e, = error due to discretization at element i f ( t )= residence time density function f ( a )= fraction of molecules at CFR outlet exiting at point a of the CFR = fraction of molecules at CFR exit leaving at point air (discretized f ) f, = linear combinator of all concentrations from PFR section of the recycle reactor fmOdel@) = weight for concentration obtained from solution of model at pth iteration F,, = mass flow rate between points ( i j ) and ( i j + l ) in the reactor segments g = reaction extent for sulfur dioxide oxidation HR,CR = set of hot and cold streams pertaining to reactor network Hp,Cp = set of hot and cold streams pertaining to flowsheet not in the reactor network H, C = set of hot and cold streams, respectively fIJ
Ind. Eng. Chem. Res., Vol. 31, No. 9,1992 2163
i, j , k = index for finite element (i) and collocation points o’, k), respectively J = objective function for reactor synthesis J, = objective function attained by a recycle reactor extension = objective value after solution of P3 $3 = objective after pth iteration of algorithm in Figure 5 1 = lower bound on Xeit N,NE = number of finite elements in the discretization procedure N p = total number of heat exchange streams K = number of collocation points LK(a)= Lagrange interpolation polynomial of degree k Lk’(a)= derivative of the Lagrange interpolationpolynomial q(a),qi, = fraction of feed entering at point a in the reactor and discretized q , respectively Qo = flow rate at entry to reactor network Q(a) = total flow rate at point a in the reactor QH,Qc = heat loads to be supplied and removed respectively by utilities R(X)= rate vector (dedimensionalized to 8-l) Re = recycle ratio t = residence time T(a),Ti = temperature profile and discretization,respectively TP = candidate pinch temperature Thin, = inlet and outlet temperatures for hot stream h tein,tYt = inlet and outlet temperatures for cold stream c T, = temperature profile within plug flow section of recycle reactor ATm = minimum approach temperature for heat transfer u = upper bound on Xerit X , X,, X , = dimensionless molar concentration vector in the CFR and discretization Xiend= concentration at the end of element i X,,X, = dimensionless concentration for components a, b X, = concentrationsavailable from solution of P3 (parametric in f ) Xo= dimensionless concentration vector (DCV) at reactor entry X, = DCV within the segregated environment Xelit= DCV at reactor exit Xu ate = DCV for fresh feed to recycle reactor extension in
&
Here, we show that (7)-(10) in P3 provide a valid approximation to the cross flow reactor. For the schematic shown in Figure 16 which correspondsto Figure 4, we have
- X a Q a = R(Xa+dAaQa+Aa + QoXoq*(a)
Xa+AaQa+Aa
(A4 Qa+Aa
= Qa + Qoq*(a)
(A-3)
where q*(a)QO is the total flow entering at point a. Also, the set (A-2) and (A-3) can also be represented by the following set of equations: Xa+AaQa+Aa
- XLQa+Aa
R(Xa+Aa)AaQrr+Aa
(A-4)
Note that (@-(lo) correspond to a particular choice for (A-4),while (A-5) is exactly (11) (recall 4 = q*(a)Q0/Qd3. It is therefore adequate to show that (A-2) and (A-3) are valid approximations to the cross flow reactor, for sufficiently small Aa. From (A-2) and (A-3),
Taking the limit as A a
-
0, we get
U / d a = R(X)+ q(a)(Q0/Q(a))(xo - Xa) (A41 where in the limit q(a) = q*(a)/Aa. Equations A-6 and A-1 are identical, and therefore the discretization for sufficiently small Aa is a valid approximation. Expression for Residence Time. The CFR residence is defined as
u = set of variables in reactor targeting formulation w,, wh = heat capacity flow rate for cold stream c and hot
7
= J m Q ( a )da/Qo
stream h, respectively y = set of variables in the reactor and energy network equations zHp = difference between heat content of sources and sinks
above the pinch point TP Greek Letters a,aij = time along length of reactor and discretization, respectively Aai = length of finite element = ai+l - ai t$i = ratio of feed flow to element i to the bulk flow entering element i Q = profit function for flowsheet design o = set of flowsheet parameters, excluding reactor-energy network (y) A = convex combinator for cross flow reactor 7 = mean residence time T R = residence time in recycle reactor
Appendix Discretization Equivalence for the Cross Flow Reactor Model. The cross flow reactor corresponding to (1) and (2) in P2 is described by the following set of equations: X ( 0 ) = xo
Also, Q i j = Qij-l + qijQo- f i j Q o , where according to the discretization in Figure 4, qij = 0 for j # 0. Thus, in the discretized model, TQO = CQijAaij = ij Qio(aii
- a101+ Qii(a12 - oil) + Q12(a13
- a121 +
for 0, 780 = (810- Q11)all+ (811- 8 1 2 1 ~ ~ 1+2 Substituting for the coefficients of the sits, we get a10 =
7
=
e..
CVij - qij)aij ij
which is (13) in the discretized model P3.
Literature Cited Achenie, L. E. K. Optimization of Chemical Reactor Networke. Ph.D. Thesis, Carnegie Mellon University, Pittsburgh, PA, 1988. Achenie, L. E. K.; Biegler, L. T. Developing Targets for the Performance Index of a Chemical Reactor Network. I d . Eng. Chem. Res. 1988,27, 1811. Andrecovich, M. J.; Westerberg, A. W. An MILP Formulation for Heat Integrated Distillation Sequence Synthesis. MChE J. 1985, 31, 363.
Znd. Eng. Chem. Res. 1992,31,2164-2171
2164
Aris,R. Studies in Optimization. I: The Optimum Design of Adiabatic Reactors with Several Beds. Chem. Eng. Sci. 1960a, 12, 243. Aris, R. Studies in Optimization. I1 Optimum Temperature Gradiente in Tubular Reactors. Chem. Eng. Sci. 1960b, 13,75. Arii, R. The Optimal Design of Chemical Reactors. A Study in Dynamic Programming; Academic Press: New York, 1961. Astarita, G. Multiple Steady States in Maximum Mixedness Reactors. Chem. Eng. Commun. 1990,93,111-124. Balakrishna, S.; Biegler, L. T. A Constructive Targeting Approach for the Synthesis of Isothermal Reactor Networks. Znd. Eng. Chem. Res. 1992,31,300. Bhandarkar, P. G.; Narasimhan, G. An Algorithm for Optimization of Adiabatic Reactor Sequence with Cold Shot Cooling. Znd. Eng. Chem. Process Des. Dev. 1969,8, 143. Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User’s Guide; Scientific Press: Redwood City, CA, 1988. Chitra, S. P.; Govind, R. Synthesis of Optimal Serial Reactor Structure for Homogeneous Reactions, Part I1 Nonisothermal reactors. AZChE J. 1985,31 (2), 185. Cuthrell, J. E.; Biegler, L. T. On the Optimization of Differential Algebraic Process Syetems. AZChE J. 1987,33,1267-1270. Douglas, J. M. A Hierarchical Decision Procedure for Process Synthesis. AZChE J. 1986,31,353-362. Duran, M. A; Groesmann, I. E. SimultaneousOptimization and Heat Integration of Chemical Processes. AZChE J. 1986,32,123-138. Dyeon, D. C.; Horn, F. J. M. Optimum Distributed Feed Reactors for Exothermic Reversible Reactions. J. Optirn. Theory Appl. 1967, 1, 1. Glaeser, D.; Crowe, C. M.; Jackson, R.Zwietering‘s Maximum-Mixed Reactor Model and the Existence of Multiple Steady States. Chem. Eng. Commun. 1986,40,41-48. Glasser, D.; Crowe, C.; Hildebrandt, D. A Geometric Approach to Steady Flow Reactors: The Attainable Region and Optimization in Concentration Space. Znd. Eng. Chem. Res. 1987,26, (9), 1803.
Glaseer, B.; Hildebrandt, D.; Glaeeer, D. Optimal Mixing for Em. thermic Reversible Reactions. Presented at the Annual AIChE Meeting, Chicago, IL, 1990, paper 7Og. Glavic, P.; Kravanja, Z.; Homsak, M. Heat Integration of Reactors-I. Criteria for the Placement of Reactors into the Flowsheet. Chem. Eng. Sci. 1988,43 (3),593-&38. Hildebrandt, D.; Glaeeer, D.; Crowe, C. Geometry of the Attainable Region Generated by Reaction and Mixing: with and without constrainta. Znd. Eng. Chem. Res. 1990,29,4+58. Kramers, H.; Weaterterp, K. R.Elements of Chemical Reactor Design and Operation; Academic Press: New York, 1963. Kravanja, Z.; Glavic, P. Heat Integration of Reactors-II. Total Flowsheet Integration. Chem. Eng. Sci. 1989,44 (ll), 2667-2682. Kravanja, Z.; Groeemann, I. E. Prosyn-An MINLP Synthesizer. Comput. Chem. Eng. 1990,14 (12), 1363-1378. Lee, K. Y.; Aris, R. Optimal Adiabatic Bed Reactors for Sulphur Dioxide with Cold Shot Cooling. Znd. Eng. Chem. Process Des. Deu. 1963,2,300-306. Levenspiel, 0. Chemical Reaction Engineering; Wiley: New York, 1962. Linnhoff, B.; Hindmarsh, E. The Pinch Design Method for Heat Exchanger Networks. Chem. Eng. Sci. 1983,38,745. Yee, T. F.; Grwmann, I. E.; Kravanja, Z. Simultaneousoptimization models for heat integration-I. Area and energy targeting and modeling of multi-stream exchangers. Comput. Chem. Eng. 1990a, 14,1151-1164. Yee, T.F.; Groesmann, I. E.; Kravanja, Z. Simultaneousoptimization models for heat integration-11. Heat exchanger network synthesis. Comput. Chem. Eng. 1990b, 14, 1165-1184. Zwietering, N. The Degree of Mixing in Continuous Flow Systems. Chem. Eng. Sci. 1959,11,1. Received for review February 3, 1992 Revised manuscript received May 11, 1992 Accepted June 3, 1992
Thermal Runaway Limit of Tubular Reactors, Defined at the Inflection Point of the Temperature Profile Shahid Bashir, Tibor Chovln, Bassam J. Masri, Ananya Mukherjee, Alok Pant, Sumit Sen, and P. Vijayaraghavan Chemical Engineering Department, The University of Akron, Akron, Ohio 44325-3906
Jdzsef M. Berty* Berty Reaction Engineers, Ltd., 1806 Bent Pine Hill, Fogelsville, Pennsylvania 18051-1501
The predicted maximum temperature difference between reacting fluid and wall to avoid thermal runaways can be exceeded in production reactors. This has been known for some time but the explanation has been lacking. The reason for this deviation was found in that the traditional approximation of the sensitivity criterion by AT 5 R P / E is correct for a limiting value at the inflection point but not at the “hot spot”, where it can be much higher. The exact expression for the limiting value at the inflection point is the total temperature derivative of the rate, and this is proven here mathematically. The total temperature derivative of a rate can be measured in a few, well-designed recycle reactor experiments. Results m e checked by computer simulation of tubular reactors. Matching to those predicted from CSTR or recycle reador (RR) measurements was excellent. The proposed interpretation explains why previously predicted limits could be exceeded in practice.
Background The runaway limit for tubular reactors has to be known for the rational design of safe and economical reactors for exothermic reactions. Experimental determination of the limit for runaway in a pilot-plant reactor is possible if this is made in one or a few t u k of the size used in commercial tubular reactors. This method is time consuming and very expensive. The studies needed to evaluate the criterion in recycle reactors are quick and much less expensive. Thus, these can be done for various catalysts and condi-
tions leaving only the final demonstration test for a pilot plant. Tubular reactors have to become sensitive before they become unstable. The sharp distinction between sensitivity and stability exists only in mathematical models. If the second-derivative term for heat conduction along the tube is left out of the model equation, only sensitivity can be demonstrated. In this case, after exceeding a limit, a very large change occurs in the temperature for a small change in one parameter, yet only one steady state is
0888-5885/92/2631-2164$03.00/0 0 1992 American Chemical Society
