546
I n d . E n g . Chem. Res. 1989, 28, 546-552
Anderson et al. (1972) and modified by Shiblom et al. (1983) for catalytically cracked products. The reproducibility of the octane number for each gasoline was fO.l. The method requires prefractionation of the total sample (containing gasoline, light cycle oil, and heavier materials) in order to remove the heavier-than-gasoline materials by back-flush from the GC column. The gasoline fraction is then separated by boiling point on a nonpolar column and the peak data divided into 32 groups of roughly similar components. Analysis of a large number of catalytically cracked gasoline samples with known research octane values allowed the use of multiple linear regression techniques to determine relative octane contributions of each group of components. Analysis of an unknown sample is then carried out by determining the relative contribution of each group and summing the product. It should be noted that the method is only applicable to samples similar to those for which the regression coefficients were determined, i.e., catalytically cracked gasolines, and cannot be valid for blended gasolines or those derived from coal or shale oil. The gasoline octane produced from the mixture of LCO and FF is close to a linear combination of the octanes obtained from the individual feeds. This result indicates that removing the recycle oil, especially LCO from the HOC feed stream, may decrease the octane rating of the cracked gasoline in the refinery. However, the octane effect may not be observable in a commercial unit since the inclusion of LCO in a riser reduces the overall cracking severity.
Conclusions It is obviously difficult to generate data from bench-scale tests that can be directly translated to predict the performance of a commercial catalytic cracker. Nevertheless, it is believed that the information presented in this article can help refiners choose between various recycle operations. Listed below are the major conclusions from this study: 1. Fresh feed had the best crackability and gasoline selectivity of the feeds tested. The conversion for cracking single feeds was in the order FF > SRO > HCO >> LCO.
2. The crackability of the single feeds is consistent with their aromaticities and average molecular weights. The more closely the cycle oils resembled fresh feed, the more easily they were cracked. 3. Mixing cycle oils with fresh feed gave yields close to or slightly worse than those calculated from weighted linear combinations of the yields from the cycle oils and the fresh feed. 4. Combining LCO with FF increased the octane number of the catalytically cracked gasoline.
Acknowledgment The author thanks Phillips Petroleum Company for the opportunity to perform this work and for permission to publish the results. Registry No. C, 7440-44-0.
Literature Cited Anderson, P. C.; Sharkey, J. M.; Walsh, R. P. Calculation of the Research Octane Number of Motor Gasolines from Gas Chromatographic Data and a New Approach to Motor Gasoline quality Control. J . Inst. Pet. 1972, 58(560), 83. Nace, D. M. Catalytic Cracking over Crystalline Aluminosilicates: Microreactor Study of Gas Oil Cracking. Znd. Eng. Chem. Prod. Res. Deu. 1970, 9(2), 203. Nace, D. M.; Voltz, S. E.; Weekman, V. W., Jr. Application of a Kinetic Model for Catalytic Cracking: Effects of Charge Stocks. Znd. Eng. Chem. Process Des. Dev. 1971, 10(4), 530. Pohlenz, J. B. New Development Boosts Production of Middle Distillate from FCC. Oil Gas J. 1970, 68(32), 158. Ritter, R. E. Is it Worthwhile to Recycle Light Cycle Oil on My Fluid Cat Cracker? Davison Catalagram 1982, 64, 5. Ritter, R. E. Light Cycle Oil from the FCC Unit. Presented a t the NPRA Annual Meeting, San Antonio, TX 1988; paper AM-88-57. Ritter, R. E.; Creighton, J. E. Cat Cracker LCO Yield Can Be Increased. Oil Gas J . 1984 (May 281, 71. Shiblom, C. M.; Fu, C. M.; White, B. J. Determination of Octane Number of Catalytically Cracked Liquid Product. Unpublished R&D Report 9551-83, 1983; Phillips Petroleum Co., Bartlesville, OK.
Received for review July 18, 1988 Revised manuscript received December 16, 1988 Accepted January 19, 1989
Approximate Dynamic Models for Chemical Process Systems Antonis Papadourakis,+Michael F. Doherty,* and James M. Douglas* Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003
Modern integrated chemical plants can be viewed as consisting of a number of smaller interconnected subsystems. Approximate dynamic models for both individual process units and subsystems are the necessary tools for preliminary dynamic studies of plantwide operability and control. A variable-order method is presented that can be used to develop simple dynamic models for such systems. T h e method is successful in retaining the dominant dynamic modes of the process, and examples are given to demonstrate this point. Most of the previous work on process dynamics and control has focused on the behavior of single processing units. The general philosophy has been that, if each unit is properly controlled, then the control of the complete process will be satisfactory. However, the current trend (Haggin, 1984a-c) is to look for ways of controlling groups of units in order to operate at the most profitable conditions. Present address: Engineering Technology Group, Rohm and Haas Co., P.O. Box 584, Bristol, PA 19007.
0888-5885/89/2628-0546$01.50/0
Since most processes contain both material and energy recycle loops, we must consider the effects of recycle dynamics if we are to obtain a realistic assessment of the dynamics and control of groups of interconnected units. The studies of Gilliland et al. (1964), Attir and Den (19781, Silverstein and Shinnar (1982), Denn and Lavie (1982),and Rinard and Benjamin (1982) indicate that recycle streams introduce positive feedback loops into the flow sheet. The recycle loops make the plant less stable, increase the gain, and make the effective time constant of the process significantly greater than the largest of the time constants Ci 1989 American Chemical Society
Ind. Eng. Chem. Res., Vol. 28, No. 5, 1989 547 for the individual process units. Thus, the behavior of a process unit when it is in a recycle loop is normally very different than when the unit is considered by itself. An integrated plant can be viewed as consisting of a number of smaller subsystems (Le., groups of units connected with one or two recycle loops). The task of performing initial studies of dynamic operability and control for complete plants can be accomplished by performing such studies on each subsystem separately. For the latter studies, approximate dynamic models for both individual process units and small plants are needed. Such models would aid in the conceptual design of plantwide operations and control policies. These could then be refined by using more detailed dynamic simulations.
Previous Work There do not seem to have been any previous studies on approximate dynamic models for small plants. However, there is a wealth of literature concerning such models for process units. The vast majority of the approximate models used by the chemical industry today consist of a gain, a dead time, and one or two time constants:
A variety of techniques have been developed for estimating K , d , T ~ and , r2,either from step response data by use of the reaction curve method (Chi-Tsung and Clements, 1982; Stephanopoulos, 1984) or from the full-scale linearized model of the process (Genesio and Milanese, 1976; Papadourakis, 1985). Some of the latter techniques, such as continued fraction expansion (Chen and Shieh, 1968), moment matching (Gibilaro and Lees, 1969; Lal and Mitra, 1974; Zakian, 1973), and Pade approximations, have the serious disadvantage that they may produce an unstable approximate model for a stable higher order system (Shamash, 1974, 1975; Hutton and Friedland, 1975). This can be caused by strong numerator dynamics (i.e., R H P zeros). The methods of dominant eigenvalue retention (Kim and Friedly, 19741, Routh approximations (Hutton and Friedland, 1975), impulse energy approximations (Lucas, 1985), and least-squares approximations (Desrochers and Al-Jaar, 1985) preserve the stability of the high-order model in the approximate one, but they can be applied only to rational transfer functions.
Models for Small Chemical Plants and Process Subsystems There are two general procedures for developing dynamic models for small plants. The statespace approach is to write the differential equations describing each process unit, the set of which forms a plant model. The other approach is to develop a simple model for each process unit (in the Laplace domain) and then combine them to form a plant model. We discuss each of these types of models below. State-Space Models. Most units are described by sets of nonlinear ODE’S of the form x = f(x,u)
rA,
1
~
Figure 1. A matrix for a small plant. uq(sl =
G,i s I
X5iS1
G5is1
u5isl
I
I XJS)
= u2isI
x2isl = ulisJ
xis)
Figure 2. Block diagram for a plant with one recycle loop.
form as eq 3. For the plant model, the A matrix is shown in Figure 1, where the Aii’s are square submatrices representing each of the process units, the lower half submatrices S i j represent the “in series” connection of unit i with unit j , and the upper half submatrices Rij represent the recycle connection (material recycle or heat integration) between units j and i. Many of the R and S submatrices will be equal to zero, denoting that no direct connection exists between the two corresponding units. If all of the Rij terms are equal to zero, Le., there are no recycle connections, then the spectrum of the eigenvalues for the p\ant will be the same as the collection of eigenvalues for all the process units in isolation. However, when recycle streams are present, the eigenvalues of the plant depend on the interconnections. Since the design of even a small plant requires the solution of algebraic equations in the order of thousands (Westerberg et al., 1979), we expect that a considerable amount of engineering effort is required to develop the dynamic equations given by eq 3. With modern executives, such as SPEEDUP, such large state-space models can be organized and maintained. However, they are best suited to testing and fine-tuning operating and control policies for complete plants, rather than designing those policies from scratch. Laplace Domain Models of Plants: Exponential Polynomials. Consider the plant with one recycle loop shown in Figure 2. Assuming that the dynamic behavior of the various units can be described by simple models (i.e., gain, dead time, and few zeros and poles), then (by applying block diagram algebra) the transfer function for the plant will have the form shown below: ~ ( s =) Gplant(s)~ ( s ) (4) where
(2)
where x’s are the state variables and u’s are the input variables. Linearizing these equations around the optimum steady-state design conditions gives the familiar linear description 8 = A x + Bu (3) If every unit in a plant is described by equations of this type, then a dynamic model for the plant has the same
where G,,(s) is a matrix whose elements are exponential polynomials of s, and the characteristic polynomial P(s,e-’) is an exponential polynomial of s given by the numerator of det [I - n:=lGi(s)]. Thus, each element of Gplant(s)has the form (GplanJij = Nij(s,e-’) /&,d Exponential polynomials have the general form
(6)
548 Ind. Eng. Chem. Res., Vol. 28, No. 5 , 1989 ni.
m
where ahj and dk are real. The order n of P(s,e") with respect to s is expected to increase as the number of units in the plant increases. It should be emphasized that, for a given plant with recycle, all the transfer functions connecting the inputs to the outputs will have the same denominator, i.e., the same characteristic exponential polynomial. The zeros of this characteristic exponential polynomial are the eigenvalues of the plant and their number is infinite (Bellman and Cooke, 1963). Transfer functions such as those given by eq 5 are difficult to use for dynamic analysis and control studies, and it is desirable to further approximate them by simpler ones which have a finite number of poles and zeros.
A General-Purpose Reduced-Order Model in the Laplace Domain Consider now eq 6. Since P(s,e") will be present in each transfer function matrix element of a given plant, it makes sense to treat the numerator of (GplanJjjseparately from the denominator, combining the results into an approximate transfer function Gapp(s)and having to approximate P(s,e-s) only once. For ease of notation, we will develop the argument for a scalar system G ( s ) = N(s,e-s)/P(s,e-s) (8) it being understood that the resulting technique can be applied to every element of eq 6. Equation 8 can be written as N(s,e-s)
G(s)=
l/P(s,eT
___ =
P(s,e-s)
(9)
l/N(s,e-s)
We approximate the two terms l/P(s,e-"), and l/N(s,e") separately, e-dNS 1 e-+ 1 z(10) P(s,e-s) Pn(s) N(s,e-S)
--- -
and combine the results to obtain the approximation
where Papp(s) = P,(s)edp. We find the unknown parameters, bo to b, and dp, and then repeat the procedure for l/N(s,e-S)and combine the results to form eq 11. In what follows, we will write dp simply as d. For stable systems, a common method for matching the transfer functions, eq 14 and 15, is the method of moments (Lal and Mitra, 1974). For unstable systems, the moments do not exist, but exactly the same formulation can be used. That is, if we expand H ( s ) in a Taylor series around s = 0 (the steady state), we obtain
H ( s ) = co
+ CIS + c2s2 + ...
where
H(')(s)
c 1. = -
i!
s=o
diH(s) H(')(s)= ds'
i = 0, 1, 2, ... (17)
Expanding Happ(s)in a Taylor series around s = 0 and equating the coefficients of the corresponding powers of s for both expansions, we obtain a system of equations relating the parameters of Happ(s) to the parameters of H ( s ) , which are given. We must note that we need to equate the first n + 2 coefficients of both expansions since Happ(s)contains n + 2 parameters (bo,bl, ..., b,, d). It is easy to show that equating the coefficients of both expansions corresponds to equating the corresponding derivatives of H ( s ) and Happ(s), evaluated at s = 0. A sufficient condition for the latter equality is the corresponding derivatives of P(s,e-s) and Papp(s) evaluated at s = 0 to be equal. Calculation of the Derivatives of P ( s p ) ) .Exponential polynomials have the general form given by eq 7 . Without any loss of generality, we require that dl = 0 and dk # 0, k = 2, 3, ..., m. Then, the derivatives of P(s,e-s) evaluated at s = 0 are given by the following equation:
Calculation of the Parameters of Papp(s), Recall that Papp(s) = P,(s)eds, where P,(s) = bo + bls + ... + b,sn. The above equation can be written as Pn(s) = E(s)Papp(s)
(11)
(16)
(19)
where E ( s ) = e-ds. From eq 19, we get
where d = dp - ds, and
+ als + a2s2+ ... + a,sm = bo + b,s + b2s2+ ... + bnsn
N,(s) = a,
P,(s)
(12)
(13)
with m < n. We justify the use of this approximation by viewing the exponential polynomials, P ( s , P ) and N(s,e"), as power series of very high order. Equation 11 is a generalized approximate transfer function which we will use throughout the rest of this work. All the transfer functions that can be obtained by using any of the model-reduction methods cited previously are special cases of eq 11. Having established the form of the approximation, we must now determine the unknown coefficients, a. to a,, bo to b,, and d. Method of Approximation. In order to obtain the unknown coefficients in the approximate transfer function Gapp(s),we write H ( s ) = l/P(s,e-S) and we approximate H ( s ) by
Setting s = 0 in eq 20, we obtain Y
E(l)(0)Papp(q-i)(O)
Pn(Y' ( 0 ) = q!C i=o
i ! ( q - i)!
q = O , 1 , 2 ,..., n + l (21)
Equating the corresponding derivatives of PaPp(s) and P(s,e") for s = 0 enables us to replace the unknown terms PapP@')(O) in eq 21 by the known terms 13(q-i)(0). It is easy to show that E(')(O)= (-1)ldL
P,(i)(0)= i!b,
i = 0, 1, 2 , ..., n + l (22)
Combining eq 21 with eq 22 yields
(14)
Equation 23 is a system of n
+ 2 equations, linear with
Ind. Eng. Chem. Res., Vol. 28, No. 5 , 1989 549 Table 111. Algorithmic Description of the Method Approximate G(s) = N(s,e-')/P(s,e-') by G,,,(s) = N,(s)e-da/P.(s),as Follows! Step 1 First, approximate the'numerator dynamics l/N(s,e-') = e-dNs/N,(s). The minimum order of approximation must be determined so that N,(s) will include the dominant RHP zeros (if any) of G(s) (see the discussion on retention of instability). Step 2 Calculate the zeroth and the m + 1 first derivatives of N(s,e-') evaluated at B = 0. Step 3 Use eq 24 to calculate dN. Step 4 Use eq 23 for q = 0,1,2,...,m to calculate the coefficients of " ( s ) . Step 5 Calculate the zeros of N,(s). If there are zeros located very close to the imaginary axis, then check whether this is an artifact of the approximation by increasing the order of approximation by one and repeating the procedure. Step 6 Approximate the denominator dynamics l/P(s,CS)= e-dw/P,(s). Step 7 Repeat steps 2-5 for P(s,e"). Step 8 Combine the results to G,,,(s) = (Nm(s)e-(dD-dN)S)/P,(s).
Table 11. Equation 23 for n 5 5
since by definition b,+, = 0. In eq 24, the only unknown is d. Solving eq 24, we can calculate d (pick the smallest positive root of eq 24 for d ) , and then we use eq 23 to calculate the b,'s in a straightforward way. The form of eq 23 and 24 for n 5 5 is given in Tables I and 11, respectively. An algorithmic description of our method is given in Table 111. In the following paragraphs, we present the properties of the method and some discussion on the retention of stability and instability. Properties of the Method. The properties of the method are listed below: 1. The method is general and can be used for the approximation of both rational and transcendental transfer functions. 2. The method retains the stability of the original transfer function. 3. The approximation preserves the dominant dynamic modes (gain, exactly; poles and zeros, approximately) of the original transfer function. 4. The approximate transfer function exhibits the same "equivalent dead time" (Matsubara, 1965) as the original one. 5. When the order of the approximation is set equal to the order of the original (rational) transfer function, the full-order model is recovered exactly. Retention of Instability. Since our method treats both the numerator and denominator dynamics of a given transfer function as denominator dynamics (see eq 91, we must be able to retain instability, in order to preserve any positive poles and/or RHP zeros. First, we consider the case where the polynomial under consideration is a regular polynomial, i.e., P(s,e-s)= Pk(s) where
Pk(s)= a.
+ als + a2s2+ ... + aksk
(25)
and Pk(s)has at least one root with positive real part, and we approximate Pk(s) = Pm(s)eds,where m
