Simultaneous modular dynamic simulation: Application to interconnected distillation columns

Simultaneous Modular Dynamic. Simulation: Application to Interconnected Distillation Columns F. S, Laganier",

J-M, Le Lann",

X. Joulia",

B. Koehret

LEAP / Analyse Fonci.ionelle des Precedes ENSIGC (INPT) • Chemin de la Loge - 31078 Toulouse Cedex, France

and

M, Morari Chemical Engineering, 210-'11 California Institute of Technology, Pasadena, CA 91125. USA

ABSTRACT Key issues related to the use of the Sim'ultaneo'us Modulal' ApP1'Oach as a basis for an efficient dynamic simulato1' are pl·esented. The lJ1'oposed stl'ategy is al,tic'Ulated 1'o'ughly at two calculation levels. First, the mod'Ule level, where the model and local control loop equations ere sol'ved all together via a global approach using an adapted Gear's 'integmt'ion code [or DAE systems. Features and problems (consiltent initialization, index p1'Oblem, discontin'uities, ...) enco'unte1'ed at this level are discussed in the proposed context. Second, at the flowsheet level. recgcle and control eq'uations al'e solved simultaneously using a DAE integrat01' along 'with Dynam'ic Sensiti'lIity Matl"ix techniq'ues 'Used to p1'Omote the convergence. Testing of this stmtegy - with compal'ison to other approaches - 'was pl'eviously performed on ftowsheet, with inte1'connected flash 'Units. It is l~xtended liere to fio'wsheets 'with highly inte1'Connected distillation columns. F'UtUl'e pe1'specti'ves and evolutions em; then discussed. KEYWORDS Dynamic Simulation; Simultaneous Modular Approach; Numerical Methods; Differential-Algebraic Methods; Sensitivity Analysis 1

Process Simulation: the Direct and Iterative Approaches

The modeling of the dynamic behavior of chemical plants yields systems of differential and algebraic equations (DAE's) that must be solved in rhe must efficient - regarding time requirements and accuracy - possible way. Like its steady-state counterpart, dynamic simulation can be tackled using two different strategies. The direct methods consist in solving the set. of equations formed by the gathering of the equations modeling all process units into one large system. In most cases, the integration of the system is reduced to the successive solutions of linear systems with a Gaussian method. The iterative methods consist in solving subsystems derived from the aforciuentionod system independently, The coupling between these subsystems is then taken into account by way of iterations on the variables appearing in at least two different subsystems. Although process simulators are usually classified as either simultaneous (equation-oriented) or modular, it has become increasingly clear in recent. years [Biegler, 1983] that both classes actually involve at least a small fraction of iterative solution, It. is true, however. that the original simultaneous approach, as used in DASP [Fletcher and Ogbonda, 1988], is basically a direct method. Nevertheless, in many equationoriented simulators, some fonu of decoupling is introduced so as t.o take advantage of specific features of the system. On the other hand, tho modular approach, for which modules corresponding to individual process units are solved sequentially, is definitely an iterative method, But in many equation-oriented simulators, some form of decoupling is introduced so as to take advantage of specific features of the system. Historically, beginning with the pioneering work of Franks [Franks, 1972], pure modular strategies were first used for dynamic simulation, However the recent trend has been towards the equationoriented approach, first used by Perkins and Sargent [Perkins and Sargent, 1982]. This strategy - or one of its variations - has also been implemented in dynamic simulators or tools such as SPEEDUP [Pantelides, 1988J, QUASILIN [Smith and Morton, 1088], DIVA [Hall et al., 1988] or ASCEND [Westerberg and Benjamin, 1985],[Piela et al., 1901], More recently, various variations of both strategies - a confirmation to the fact. that all problems are not solved yet - have been presented under different acronyms such as two-level approach [Ponton and Vasek, 1086], modular integration lCunently on leave at Caltech, Prof', MOHAll I', COlli 1'01 (lronl' 2Conespondcnce should ~J(' sent to : ph ....", 1l1-52-~)2-1,] (,,17501' 'It;U), fax 61-.53-69-86

5287

S288

European Symposium on Computer Aided Process Engineering-2

[Liu and Brosilow, 1987], parallel modular approach [Hillestad mid Hertzberg, 19861, sequential-clustered strategy [Fagley and Carnahan, 1990J, structured mixed-mode equation-oriented simulation strategy [Gani et al., 1990J or equation-oriented [Juslin et al., 1991J. However, as mentioned by Biegler [Biegler, 1983], it may be noticed that very few attempts have been made to use the simultaneous modular framework - whose suitability in steady-state applications is well known - for dynamic applications. Basically, the main arguments in favor of each approach can be summarized as follows. The modular approach is structurally appealing to the process engineer, as it mimics the behavior of physical flows of matter and information, It can use prior knowledge of individual units in an efficient manner. In particular, it is possible to use f01' each unit an integration code that is well-adapted to its usual dynamic behavior, which may thus dispense us with unnecessarily elaborate strategies. In return, some kind of synchronizing procedure has to be thought up for the management of recycled streams and the correct integration of all units in the process. Simultaneous approaches, as long as they involve no iterative strategies, do not require such a synchroniziug procedure. More elaborate schemes, however, do need it. These try to exploit differences in the response times of groups of variables (flowrates and compositions for instance) and create a partition of the original system in an attempt to decrease the stiffness inside these subsystems, and thus to allow for the usc of explicit integration formulae. On the whole, the most convincing argument in favor of either approach relates to the solution of the linear system to which we referred earlier. Direct. methods using a form of Gaussian or LU strategy to solve it have a o(n 3 ) complexity, when' /I is the size of the system. Taking advantage of this cubic dependency, iterative methods can become more efficient, since they handle much smaller subsystems, as soon as the required number of iterations can be bounded by a reasonable value. While it is true that other factors may influence the efficiency of each approach (such as the multirate integration aspect of the modular approach for instance), it is the opinion of these authors that they do not affect significantly the relative computation costs of the different approaches. The strategy we present. in this paper. called the Simultaneous Modular Approach, relies mainly on three ideas that will be developed below :

eA two level integration strategy, where modules, corresponding to individual units, are solved at a. lower level, while recycled streams and dclocalized control loops (regulatory control loops for which the controlled variable and the manipulated variable belong to two different units) are integrated at all upper level. e A variable time horizon strategy based on the control of truncation errors of the upper level variables, that aims at minimizing the number of upper level iteration while taking advantage of different latencies in the units, eA possible interaction in both directions between the upper and lower levels, so as to generate an efficient convergence operator for the flowshcet (upper) level variables using information created at the module (lower) level. Results presented in an earlier paper [Lagauier et. al., 1991], couceruing the dynamic simulation of interconnected flash and mixed tank units, had already illustrated the feasibility and possible efficiency of this approach. While not very industrially relevant, these studies had the interest of displaying sharp dynamics (dynamics that are not. found in the case of the smoother distillation column) designed to place our approach in difficult conditions. Additional results, as well as a more realistic application example, are given here. We focus this paper on a single process of two interconnected columns as it displays the basic features of most examples we studied. The first section of this paper discusses the model we used for the distillation column, and gives a hrief review of the integration method that is used to solve it. The next, section covel'S the description of the simulation package, and goes into the details of its main features: the Coordinator, the time horizon strategy, and the generation of the flowsheet convergence operator. Finally, the last section coutaius a brief description of the illustrative example we provide, along with a discussion of the results which have been obtained. We conclude with an outline of the further research that must be conducted to improve the performance of the proposed approach.

2 2.1

Dynamic Models and their Solution Dynamic Model for the Distillation Column

Multistage multicomponent separation has been studied extensively, and in recent years a fair number of papers dealing with the subject have been published. Rigorous stage by stage models are now available, and the one used ill this work relies heavily on the one recently proposed by Le Lann et al. ([Le Lann et al., 1990]). The basic assumptions that are used an' the following : eThe vapor holdup ou each tray is neglected. This is a reasonable assumption as long as the pressure in the column is not too high (ie., P < 10 ann, sec f01' instance [Choe and Luyben, 1987]). eThe pressures ill the {"01un11l remain constant throughout the simulation, This is quite accurate as long as the pressures are not 100 close to vacuum and when there is no floating pressure control. eThe liquid phase and the vapor phase an' assumed to he both perfectly homogeneous. eThe thermodynamic equilibria arc assumed to be instantaneous, i.e. the liquid on one tray is always in equilibrium with th~~ vapor leaving that tray (with the option of taking into account the tray Murphree efficiency). Moreover, the temperature of the liquid leaving a tray is the same as the one of the vapor leaving the same tray. This is acceptable (july if the changes in temperature from one

European Symposium on Computer Aided Process Engineering-2

S289

tray to the other are not too large (see [Seader, 1989] for example). The variables whose dynamics are computed me the vapor and liquid flowrates, temperatures, vapor and liquid compositions, liquid holdup and enthalpy on each tray. They are computed from the following equations: eMaterial Balances for each component. e Heat Balance. eTotal Material Balance. e Normalization Equation for the compositions, eEquilibrium Equation for each component. eLiquid Holdup Model (constant volume, tank model or Francis weir formula). e Liquid Enthalpy Model. In addition, regulatory control loops can be added on any temperature or molar fraction in the column by acting on one or two of the operating variables (heat dut.ies, reflux or reflux ratio, distillate f1owrate). The steady-state simulator PROSII'I'f [Koehret 1987] is interfaced with our programs so as to provide the initial steady state from which our dynamic simulations are started. The sets of equations we alluded to above can be written in a general form as : F(t,'y(t),y(t),'u(t))

=0

(1)

where F : ~ x ~" x ~.. x ~" - t ~" is the set of equations describing the model, y E lRn the vector of variables and u E 1RP a vector of inputs. The general features of these systems in chemical applications are now pretty well known, and the main ones are the following: e Moderately to highly stiff systems. The time constants of the variables involved in the system spread over many orders of magnitude. As a result, explicit. numerical methods perform very poorly on these systems, as their time steps tend to become unacceptably small. e High index systems. The index of a system of DAE's is defined as the number of times the system or part of it must be differentiated so that the time derivatives of all variables can be expressed explicitly as a function of the variables. It is known that while usual numerical methods behave fairly well on index 0 or 1 systems, they tend to either give wrong results or fail on higher index problems. e Highly non-linear equations. The equations modeling thermodynamics and hydrodynamic correlations are highly non-linear, which increases the coupling between variables. e Large but sparse systems. As was already mentioned, the system of equations modeling a process can reach very large sizes, especially if distillation or absorption columns are involved. However, t.hose systems are very sparse, as only a few of the variables appear in each equation. e Discontinuities. These appear frequently iu dynamic systems, as the number of phases may change when temperature, pressure and composition vary for instance. However, this aspect of the problem won't. be tackled in this paper, as we consider ouly two phases system. However, it should be noted that problems with discontinuities are very important for start-up and shut-down procedures, for cyclic mode operation and when dealing with batch, semi-continuous or discrete systems. 2.2

Solution of Systems of Differential and Algebraic Equations

While other met.hods are known and have been studied intensively ( semi-implicit Runge Kutta methods, for instance [Prokopakis and Seider, 1981]), codes based on the Backward Differentiation Formula of Gear [Gear, 1971] such as LSOm [Hindinarsh. 1974] or DASSL [Petzold, 1983] are more popular. It should however be kept in mind that these solvers were not developed for general DAE systems, and as such should not be expected to give systematically accurate results for any arbitrary problem. Moreover, index 2 problems do not belong to the class of problems that are expected to be solved by these codes (see for example [Bronan et al., 1989] for a discussion of the index problem). Therefore, special care is recommended when using these codes to solve chemical engineering problems. Given data on the past behavior of tho solution. the system at time t is solved via a predictor-corrector scheme, where the derivatives of the variables are extrapolated from the values of these variables by : Ji(t) = Cj (y(t) -t/J) (2) Here, Cj is a coefficient that varies with the order of the method and the integration time-step, and'" is a function of past values of y . (2) is then substituted into (1) to yield a set of algebraic equations that are solved iteratively: F(t,y(t),Cj(y(t)-rP),u(t))

=0

(3)

The predicted value y(O) and it.s time derivative !itO) are obtained by polynomial extrapolation based on previous timepoints, The corrector equation, for a fixed timestep, is solved using a modified NewtonRaphson iterative scheme, given by y(m+l)

111

= y(m l _

(l

J- 1F(t. ylm), (j(Y("') -

= 0.1, 2, . _.

~',

II)

(4)

where J is the itemtio'fL matl"i.x, J

DF =-of +Cj-o. oy y

(5)

S290

European Symposium on Computer Aided Process Engineering-2

and a is a scalar constant chosen to speed up the rate of convergence of the iteration when J is the iteration matrix from some previous timepoiut. For it more detailed formulation see [Brenan et al., 1989]. The stepsize and order selection strategies in DASSL tend to favor sequences of constant stepsize and order. There is an initial phase where the order and stepsize are increased at every step, and after that the order is lowered or raised depending on whether the estimates for the leading error term in the remainder of the Taylor series expansion, IllIk-Iylk-llll, IIhkylkJII, Ilhk+l y(k+l)lI, and IIhk+ 2 y(1c+2) II form an increasing or decreasing sequence, respectively. The new stepsizc is chosen so that the error at the new order k, estimated as if the last k + 1 steps were taken at the new stepsize, satisfies the condition MIlY,,+1 - y~.oJllI

:s:

1.0

(6)

where 111 is a function of the order and the step length of the previous steps which bounds the error estimate taking into account the interpolation error and the local truncation error. To control the local error the current stepsize is rejected whenever the condition (6) is not satisfied. Note that condition (6) requires that the estimated integration error be less than the requested integration error tolerances, as reflected in the weighted norm of the predictor-corrector difference

11.//11

=

1 ---;

N

'"

i\' L

..

;=1

Yi 2 (uit, -)

IIY,,+! - y~~ll1, where lIyll is given by (7)

'

where N is the system dimension. and 11't is the weight vector rofiectiug the relative and absolute error tolerances : (8) wt[i] atol[il + l·tol[i] * y[i] where atol[i] is the absolute errol' and I·tol[i] the relative error. In addition, in order to exploit the sparse structure of the system, the sparse linear solver Sparse1.3 [Kundert and Vincentelli, 19881 has been interfaced with DASSL. It was rationalized that the banded option available in DASSL WOUldn't take full advantage of the sparsity of the system, and the results we obtained do indeed justify this.

=

3

3.1

The Simultaneous Modular Approach The Two-Level Strategy

The strategy which has been developed for the SMA rests on the realization that in the context of modular approaches, the dynamic behavior of the torn recycle streams is governed by an implicit (since hidden in the models of the different units) DAE's system, which can as such be solved with an adapted integration code. Moreover, it is based on the assumption that. an efficient error control at each mesh point of the interval is enough to ensure consistency at any other point of that interval. As we mentioned earlier, each unit in the process is modeled as a module that can be solved, given its input, for its state variables. The modules are solved sequentially over a time interval, the output of one unit being used as the input for the next unit. Since we consider dynamic units, the value of the output will in most cases be different at each point of the interval, and an appropriate representation must be used to transfer accurately these time functions from one unit to the other, This is done by way of Lagrange polynomials that interpolate the functions at P of their previous computed values, Moreover, if a stream leaving a unit is recycled as an input to an earlier (in the sense of the fiow direction) unit, we have to "tear" this stream and guess its value so as to be able to compute tho Iat.tcr unit. The value computed as an output to the former unit is then compared to the guessed value, and iterations are performed until convergence is reached. Now, since the computed value is deduced from the solution of the units that are affected by the torn stream, it follows that there is an (unknown) system of DAE's that the torn stream must satisfy. It is clear then that the equations for the tom stream can be solved with a predictor-corrector method, where the corrector residuals arc evaluated at. each iteration by integrating the whole process on the time interval defined by the time step used for the integration of the torn stream. So, we end up with a two-level integration strategy, where the modules are solved over a given time horizon with as many steps as necessary at a lower level, while all torn streams and delocalized control loops are solved simultaneously at an upper level. The integration of all modules is synchronized via what we call after Liu and Brosilow [Liu and Brosilow, 1987] the Coordinator, a modified integration code derived from DASSL, which handles the solution of the upper level variables. 3.2

Error Control and Time Horizon

Since the upper level variables are solved with an integration code, two enor controls are performed on these. The convergence error checks ensure that the convergence at the mesh points for the upper level variables is achieved. While it is deal', in the context of steady-state modular simulators, that the convergence check should be uiore stringent in the modules than it is for the flowsheet variables, it doesn't follow that it is necessarily tilt' case for the dynamic simulation. On the opposite, our results suggest that it is likely that the same error tolerance should be used at. both levels. The convergence between two mesh points is ensured by keeping the truncation (local) error on the flowsheet variables bounded. This truncation error serves as a basis for the computation of the time horizon, the time step for the fiowsheet variables, which is derived so as to minimize the number of flowsheet. iterations. It turns out that the parameters used in DASSL for the modules are fairly reliable for use at the flowshcet level. Expcrimeuts in rhaugiug these parameters did not improve the results on

European Symposium on Computer Aided Process Engineering-2

S291

the average. It should be understood that the time horizon length is governed mainly by the dynamics of the tom streams, and the ability of the predictor to give a good guess for the initial value of these streams over the next horizon. It follows that some units may require many steps to integrate over a given horizon, while some other may be integrated with only a few steps. Also, since the the horizon is computed using the truncation errol', it is clear that the order p used for the Lagrange interpolation polynomials should be equal to the order of the integration formula at the upper level. 3.3 3.3.1

Flowsheet Convergence The Flowsheet Convergence Operator

To solve the corrector equation at the fiowsheet level, we have seen earlier that a Jacobian matrix must be provided, 01' at least an approximation thereof. Also, we have seen that the tom streams that are solved at that level are governed by a DEA's system that is staticized by the corrector formula as : Z - G(Z)

=0

(9)

=

I - ~,where ~ measures the influence of a change of the torn stream variables on the torn stream, or to put it

where Z is the vector of torn stream variables. The Jacobian matrix is then J

another way, ~~ is the sensitivity of the torn stream to itself. Now, if we realize that the sensitivity of the torn stream to itself can be decomposed as a combination of the sensitivities of output streams of the modules to their input streams, we can generate very easily the flowsheet convergence opera.tor by using the Dynamic Sensitivity Matrices (DSM) of the modules. As an illustration, consider the process depicted on figure 1. If Cl is chosen as the torn stream, we want to solve for each mesh point the equation Cl - G( Cl), where G is a function of the three unit models of the process, The sensitivity may be computed by analysing how a perrurbariou on Cl propagates back to itself, namely paths Cl- C2 - C3 and Cl - C5 - C3 - C1. Thus the sensitivity of torn stream Cl to itself can be expressed as a function of the modules DSM :

ttl

= SB(I,2)·S,\(2,1) + SB(I,5).Sc(5,3),SA(3,1)

. S(I,I)

(10)

where Sx(i,j) is the DSM of stream Ci to stream Cj in module X. CI

C2

co

C5

C3

C4

Figure 1: Computation of the Flowshect Operator Since these module DSM are gencratrd, as w« will see in the next paragraph, at the module level, we thus use module level information to promote convergence at the Howshcet level. 3.3.2

Dynamic Sensitivity Matrix of a Module

If the model for module i is : F;( t . U;( t) , u"; ( t) .u; ( t))

=0

(11)

then the DSM of this module is found by computing the following matrix product: S;(y;, lI;)

u, = -J; .-/-1

(12)

( lIj

where J, is the Jacobian matrix of module i computed for a time step equal to the horizon length, and ~ is the derivative of the model equations for module i w.r.t, the input variables lIi. Note that the DSM generated in this way contains the sensitivities of all variables of the module to all its input, and not the stream-stream sensitivities only. Since the module Jacobian matrix has usually beeu last evaluated for a time step different from the horizon + Cp ~,f' and LU dccompositiou must be performed again (though all length, the computation of u. ««, P individual terms y, !ht and C1' are already computed). J

fl

fl,

CACE17 suppl-T

tYj

S292

3.4

European Symposium on Computer Aided Process Engineering-2

Alternate Choices for the Convergence Operator

While it has been shown that an efficient convergence operator may be computed at very low additional cost, it may be argued that a simpler operator may be used for the flowsheet level, especially if we keep in mind the fact that the initial guess deduced from the earlier values of the variables will usually give a very good estimate of the real value. To a certain extent, the identity matrix (leading to a successive substitution method) may be used with pretty good results. It has actually been seen on the examples we studied that the results obtained with the identity matrix are usually so good that the use of the DSM is not really required, However, it may tum out that they do improve convergence significantly in other applications, especially with stiffer or 1I1On~ non-linear systems, or in the case of sharper transients. They may also be particularly helpful whou dealing with delocalized control loops or flowsheet level constraints. In the example we present, DSIVI's are computed numerically, which may admittedly be an inefficient and unreliable way, In future works, it should be obtained analytically, so as improve the overall scheme efficiency. The systematic use of the DSM concept in the context of dynamic simulation at a flowsheeting level is relatively new [Laganier et al., 19911' though it has already been used successfully in steady-state simulation for optimization purposes [Wol xert et al., 1991], for parameter estimation in DAE systems [Caracotsios and Stewart, 1985]' \Dillarddlo et al., 1992], and in the context of optimization for DAE systems with optimal control app icatious in mind [Fahrat et al., 1990bJ, [Fahrat et al., 1990aJ. The key point remains in the way the DSM is obtained and how it is exploited efficiently. 4 4.1

Application: Interconnected Distillation Columns Description of the process

For the purpose of this paper, W(~ studied a process with two interconnected columns, with material side-streams on each column, originally proposed by Hess and Holland [Hess and Holland, 1976]. The feed stream is a mixture of 7 hydrocarbons, upon which a disturbance (in flowrate, temperature or composition) is applied. The form of the disturbance can be either a valve opening, a step change or a sinusoidal perturbation. Both column have 52 trays, and the nominal feed flowrate is 100 mol/hr. We call column 1 the column on the left of figure 2 and column2 the other one. Other data may be found in the referenced paper. CI

10

FEED

51

21

45

C3

21

52 31

C2

C4

Figure 2: Interconnected Columns The thermodynamic model is interfaced from PROSIJ\·!. The hydrodynamic models are the following: Francis formula for the trays, constant volumes for the rehoiler and condenser tank. The results presented here conespond to a 10% ramp decrease in the feed flowrate, For those results dealing with the tolerance ratio between fiowsheet and module levels (see below), we used a similar but smaller example (32 and 10 trays respectively] and studied the dynamics created by a 10% decrease of the reflux ratio of colunml. 4.2

Results Fitting

In order to check the accuracy of the results, since no dynamic data are available, the results have been compared with those obtained by solving the whole system of equations with DASSL with a very tight tolerance, These are displayed on the plots as dashed lines to which the computed points should be compared, and are designated henceforth as the reference. As it is clear on figures 3 and 4, the agreement between our results and the reference can be made as good as desired by manipulating the error tolerances. Stream C2 displays inverse responses, while stream C3, which is reached by a damped perr.urbnr.ion iuuch later. has a monotonic response. 4.3

Influence of the Choice of the Torn Stream

Traditionally, tom streams are selected with tools from the graph theory, that is, they are chosen by considering the network conuectious of the prOCt'SS, so as to minimize the number of such torn streams.

S293

European Symposium on Computer Aided Process Engineering-2

6

li.O·3

--- Reference

0

--- Reference

Computed points

0

Computed points .~.,

.J"'!t-

.....- . . -~J. . . . . .

5

0.04

0:

0

.~

I '(;

'ij :>

~ 0.03 i5.

/0

"._~..

9~

t

~

2

••.•. C3

/

I

l;j

0

-1 0

,..

••' , ••0'

~0'01~ >

1000 2000

3000

4000

5000 6000 7000

o

8000 .0.010

~.

,,'

I

0.02

~

~

1'''

..,....

Il,r.>"'

3

V,I

.!

Q

::l:

4

ol::

..!to.....,.-.0.9-....

~

•.••

0-'0" Q4"""'''

.,.J'

1000 20
Time (s)

Time (s)

Figure 3: Results fitting for tol=le-3 and torn strcam=C2 -.. Reference

0

... Reference

Computed points

0

Computedpoints

5 4

C3

3 2

1000 2000

3000 4000 5000 6000 7000 Time(s)

8000

1000 2000

3000

4000 5000 6000 Time(s)

00

Figure 4: Results firring for tol=1L'-3 and torn strcam=C3 In the context of dynamic simulation, other crit.oria should also be taken into account, such as the time constants of the units and the effect input disturbances have on their output, In the case of Oul' example, either of the two connecting streams can be torn. However, it seems advantageous to teal' stream C2 if one is interested in the accuracy of the flowrate responses. Indeed, these have very fast responses (time constant is around 400 s), and thus will dictate the size of the horizon during the first part of the simulation. Given that the columns in am example are parametrized with the distillate flowrate (perfect control), the flowrat.e of stream C3 remains constant throughout the simulation. If C3 is used as torn st.ream, t.he horizon size is decided according to the evolution of the compositions in C3, and thus increases rapidly at the beginning of the simulation, which implies less accuracy for the rapidly changing flowrates in the columns. Moreover, if C2 is chosen as the torn stream, the timing results are highly dependent on the initial horizon size which must ]w chosen more accurately than in the case where C3 is torn, for which the dynamics appeal' only after a significant time lag. 4.4

Horizon Changes and Multirate Integration

It can be seen in figure 5 (a) that the evolution of the horizon size is highly dependent on the choice of the torn stream. As we mentioned previously, if C3 is tom, then the horizon size is governed mainly by the composition dynamics. Moreover, the disturbance being applied on column l , there is a time lag during which stream C3 is hardly affected and during which the horizon increases steadily. On the other hand, if C2 is picked as the tom stream then the disturbance affects it almost immediately, which accounts for the fact that the horizon remains quite small for a while and then increases rapidly. The stepwise evolution of the horizon size is due to the specific strategy used in DASSL for the stepsize selection, which was kept unchanged in the Coordinator. On the plot of figure 5 (b), it. becomes apparent that. the rnultirate effect does indeed exist, though it quickly disappears as the simulation procoixls. Also, it must be noted that. it is in close dependence to the ratio between the tolerunccs of rlu- fiowshevt and module lovel, and is actually quite insignificant if this ratio gets very dose to 1. The results presented on the plot. were obtained for a 3/1 ratio, which we show next is not a very appropriate rat io.

4.5

Tolerance Ratio, Accuracy and Timing Results

It is commonly accepted that the ratio between the fiowsheet error tolerance and the module error tolerance should be larger than 1. If one does not give it. more thought, it seems even likely that the

S294

European Symposium on Computer Aided Process Engineering-2 (8) - Slrearn C2 tom • • • Stream C3 torn

900 800

(b)-Colwnnl ·- - Column2

I:;

700

....

/1

600

':::'500

S400

'il

/

::r:

300

200 100 /

00

.:

!

""! i

i ("--

.

if 2000 3000 4000 50'00

6(iOO-1000-11oo0

00

SOO.----d.".-.--J;;~~;;r=:;r;;;;;r=:;;:;t;;;;r=:;r;;;;;:::::;:;!~

Time (s)

Time (s)

Figure 5: Influence of the tom str eam on the horizon ami average number of time steps in each module per time horizon method will fail if a ratio lower th an 1 is used, However, it nmy also be argued, from the viewpoint of the modules, that this ratio should be 10Wl'1' than 1, since the output of a module can be computed with a given accuracy only if its input is provided with at least the same accuracy. Thus a ratio equal to 1 may actually be optimum . For lack of time, the studies displayed au figures 7 and 6 were conducted on the smaller example mentioned earlier. On these, it. app ears tha t no significant gains in CPU time can be obtained by increasing the ratio over 1 (JUol designates the tolerance used at the flowsheet level). However, the accuracy of the results does change dramatically, with the ratio. The plots we present show the evolution of the temperature of the connecting stream which is not used as torn st ream, compu ted as an output of its relevant module, so that these plot s reflect th e iufiueuce of both tolerances. Torn Stream : C2 . •, ~ .. .. " '.' , . ..'• .., .

900 :.- - .- -,. ~ -

lOoor-~- '~ ~~- -' -~ -" "

950

~lou ulc

lISO:

Tolerance : 5c· 3

900

lIoo

850 \"

S800

~

..

750

600

37S0

\ ...,

..

'\

~ 700 u 650

Tom Stream: C3 ..~~-~~~~~-~..,....~~."

700

6S0f

:::>

\ ...

Module Tolerance : le·2

~

eJ

Module Tolerance : 1e·2

'\

..

......

.........

550

5q
:q':l_- -~ ,~=::;:::==;::;::====~ u·>

lO·'

10'

Tolerance Ratio (F1owsheel/Module)

Figure 6: Influence of the ratio between fiowshect and module tolerances on CPU More surprisingly, it. app eal's that the accuracy of the results can be increased by reducing the ratio even further, past the value L In this case, the tigh ter tolerance on the torn stream creates an increase in the number of fiowshect mesh point s, which ill turn increases th e CPU time because of the additional stop points it imposes at the module lev el, However, whatever the value of till' ratio. the simulation with the S!\lA of these particular applications were about three times slower than with a direct approach. Thi s can be explained by the two following facts: .At least 2 flowsheet iterations per horizon are required to improve the stability of the method. This means that the whole problem is actua lly solved at least twice (not much more, really, as it turns out that the average number of Howsh ect iterations lit's between 2 and 3). • The integration of the modules is forcingly stopped a t each mesh point of the flowsheet integration, again for stability reasons. This accounts for III ore time steps than may be necessary for the global approach. In any case, it is our belief that gains in CPU time Call be achieved only if more than 2 units are involved in the example. If we accept the assumption that most of the computation time is spent in the solution of the linear system, then the total CPU time can be expressed as a function of the average number of flowsheet iterations 1\Ti 1c n the number of modules l· and the total size of the system N :

European Symposium on Computer Aided Process Engineering-2

S295

383

1000

Torn Stream: C2

Tom Stream: C3

Module Tolerance le-Z

Module Tolerance 1e-2

2000

T.5M . \

= C.5· A.

i t e l, •

. N 3 k . (k)

(13)

where C s is proportional to the average number of steps taken in the modules. The CPU time for the direct approach is : TClLo

=

J

1\iJ

(14)

CS'H

The SMA will be faster than the direct approach only if

¢}

k

>

ll'iter .

C·'s

Cs

(15)

If we take the approximate values - suggested by the ones observed III the simulation we performed Niter = 2.5 and Cs = 1.5 C~, we find that k umst be at least larger than 2. 5

Conclusion and perspectives

Some key features involved in the clcvolopinout of a new suuultancous modular approach for dynamic simulation purposes have been enlightened in this contribution. The most prominent of these features is that the proposed approach is based on a two level strategy so that the characteristics of each level can be taken advantage of. At the module level, most of the mathematical procedures developed for the equation-oriented approach can be used, with the additional flexibility of introducing pertinent physical intuition at different stages of the problem. This allows for high-fidelity models and the solution of very large and complex problems. Discontinuity problems, index difficulties, time and state events detection, problems with fast transients may be solved locally without interfering at. the fiowshect level. For that reason, the possibility to diagnose and locate the type of the encountered difficulties is enhanced. Using the same type of DAE integrators (or similar ones) at both levels ensures a gl"P'H. consistency. However, this does not preclude the use of different implicit integrators for clifforcut modules - including tho coordinator - provided the integration scheme involves an iteration matrix inforuiativc enough to generate the DSIVI's of the modules. Moreover, the difficult problem of' finding a consistour init ial s\'t of variables when dealing with DAE systems is simplified. At the fiowsheet level, different. classes of problems which do not exist in the context of equation-oriented simulators have to be investigated : choice of the tom streams, time and state events management, interpolation order, numerical iiuproveiueuts in the coordinator scheme, use of a reduced-order dynamic operator, etc. These were apparently solved successfully, as the approach was able to solve both fast transient problems as in the case of intercounectcd flashes [Laganier et al., 1991] and very large scale, highly interconnected problems, as in the case of coupled distillation columns. As we mentioned earlier, the dynamic sensitivity analysis uSNI to improve the flowsheet COli vergence has been already applied to different. problems. Its usefulness to the siumltaueous modular approach should appear more dearly in the context of dynamic optimization. The approach has been tested on tlowshcors with inrcrcouuccted distillation columns, with systematic comparisons with tho global approach. Givou til<' results that wen' obtuiucd, this approach may form a base for the development. of a general purpose dynamic simulator. At the module level, introduction of new specifications in the standard module description will be iuvcst igated. Thus, very general dynamic models, including sets of options to IJ\' chosen froiu for each specific application, can possibly be created from a standard base (the existing standard modules). Methods oflines (fixed grid) which have been used to solve distributed parameter models [Heydwoiller et al., 1977], [Bacchus et al., 1992] have appeared to be suitable for the proposed approach. A future topic of research is to investigate the use of methods

European Symposium on Computer Aided Process Engineering-2

S296

of lines with an adaptive moving grid [Petzold, 1988],[Bl"Cnan et al., 1989]. At the fiowsheet level, in cases of very highly interconnected problems 01' when dealing with systems with strong coupling among compositions, flowrates, temperatures ami especially pressures, where the global approach is obviously more efficient, the involved modules can be aggregated into a block-module, on which the dynamic sensitivity analysis is performed. From a more general perspective, in connection with ideas developed previously - conception of a general purpose continuous/discontinuous simulator [Daubas et al., 1991J, simulation of multisequence batch distillation processes [Alber et al., 1991J - the application of this strategy does not need to be restricted to the dynamic simulation of continuous processes and can be extended to the cases of discrete or semicontinuous problems (for instance, start-up and shut-down procedures and batch processes). This approach makes it. possible to take into account all the important work done in the modular context together with the impressive research advances that have been made along the equation-oriented lines. Though substantial work remains to be done, pursuing this course seems a promising and exciting endeavor. 6

Acknowledgement

This work was completed as part of F. Laganier's Doctorat thesis (Institut National Poly technique de Toulouse). We gratefully acknowledge Rhono-Poulenc Industrialisation and the Rhone-Poulenc Company for their financial support. References Albet, J., Le Lann, J., Joulia, X"mHI Koehret., B. (1991). Rigourous simulation of multicomponent multisequence batch reactive distillation. In Comzrukr-(h"i.ented Process Engineering, Proceedings of COPE-91, October 14-16, Barcelona, SZla·jn, pagt's75-80. Bacchus, P., Jacob, J., and Le Lauu, J. (1992). Simulation cl'un filtre immerge de denitrification: modelisation et resolution numerique, Technical report, DEA reports, GAFGP #91-13 and #92-03. Biegler, L. (1983). Simultaneous modular simulation and optimization. In 2nd International Conference on Foundations of Comp·utel·-A·jded Process Desiqn, pages 369-409. Billardello, P., Joulia, X., Le Lanu, J.. Delmas, H., and Koohret, B. (1992). A general strategy for parameter identification in differential-algebraic systems, In Buropem, S'ymposium on Computer Aided Process Eng·jneel'ing (ESCAPE-1), 26-28 Mug, Elsinore, Denemark (SuPIJlementary Vol.), pages 165172. Brenan, K., Campbell, S., and Petzold, 1. (1989). Numerical Solutuni of Initial Value Problems in Differential Algebm'ic Equations. North-Holland, New York. Caracotsios, M. and Stewart, W. (1985). Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Compo Clietn, Eng., 9 (4):359-365. Choe, Y. S. and Luyben, W. 1. (1987). Rigorous dynamic models of distillation columns. Ind. Eng. Chem. Res., 26:2161-2164.

Daubas, B., Pingaud, H., and Koehret, B. (1991). Developpement d'un simulateur de precedes discontinus, semi-continua et continus. In Rlicents Proqre« en Genie des Procedis, 3eme Congres Genie des Precedes, Paris, France, volume 5 (13), pages 307--312. Lavoisier. Fagley, J. and Carnahan, B. (1990). The sequential-clustered method for dynamic chemical plant simulation. Compo Chem. Eng., 12 (5):161-177. Fahrat, S., Czernicki, 1\1., Pibouleau, L., and Domenech, S. (1990a). Optimization of a multiple fraction batch distillation by non-linear programming. A/ChE Journal, 36 (9):1349-1357. Fahrat, 5., Pibouleau, L., Domenech, 5., Alber, J., Joulia, X., and Czemicki, M. (1990b). Optimisation d'une operation de rectification discontinue multiconstituant multifraction. In Colloque Inf01matique et Precedes Discontinue, Paris, France, pages 11--16. Fletcher, J. and Ogbonda, J. (1988). A modular oquation-orionted approach to dynamic simulation of chemical processes. CO'Inp. Clunn, Eng .. 12(5):'101-405. Franks, R. (1972). Modding a:nd Siniulation. in Chemical Engilleel'in9. Wiley, New York. Gani, R., Perregaard, J., and Johansen, H. (1990). Simulation strategies for design and analysis of complex chemical processes. Trans, I. Cliem., 68 (A):407-41 7. Gear, C. (1971). The automatic iutegratiou of ordinary differential equations. Communications of the A.C.M., 14 (3):176-179. Hess, F. and Holland, C. D. (1976). Solve more distillation problems - part 6 : evaluate complex or interlinked columns. Hyd1'Oca.'fboll P1'Ocessing, June: 125-131. Heydweiller, J., Sincovec, R., and Fan, L. (1977). Dynamic simulation of transient systems described by distributed and lumped parameter models. Comp. Cliem. Eng., 14:813-869. Hillestad, M. and Hertzberg, T. (1986). Dynamic simulation of chemical engineering systems by the sequential modular approach. Comp, Clicm, Eng., 10 (4):377-388. Hindmarsh, A. C. (1974). LSODI : ordinary diftorential equations system solver - rep. ucid-30001, rev.3. Technical report, Lawrence Livermore Laboratory. Holl, P., Marquart, W., and Gilles, E. (1988). DIVA - A powerful t.ool for dynamic process simulation. Camp. Chern. Eng., 12 (5):421-426. Juslin, K., Kaijaluoto, S., Kalitvouzeff, D., Kilakos, A.., and Lahdenpera, E. (1991). An equationoriented software package for dynamic simulation of chemical processes. In Computer-Oriented Process Engineering, pages 405-409, L. Puigjanor, A. Espuna, ed. Elsevier. Koehret, B. (1987). COnCelJtion d"un .~im·ulat(~-Ul' de In·oddJ.~. PhD thesis, I.N.P. Toulouse.

European Symposium on Computer Aided Process Engineering-2

S297

Kundert, K. and Vinccntelli, A. (1988). A sparse linear equation solver - User's guide. Technical report, Dept. of E.E. and C.S., University of Berkeley, Califoruia. Laganier, F., Le Lann, J., Joulia, X., and Koehret , B. (1991). Use of the dynamic sensitivity matrix concept for the efficient solution of dynamic prOCl'SS simulation. In Pvoceeduiqs of The IV International Symposium on Process Systt:m Eng'ineeTing, Mont(~bello, Canada, pages 000-000. Le Lann, J., Alber, J., Joulia, X., and Koohret, D. (1990). A nmltipurpose dynamic simulation system for multicomponents distillation columns. In EFChE S1j'mlJOs'i'um ComChem '89, La Haye, Computer Applications in Chemical Engineel'ing, ll<'\gl'S 355-361. Liu, Y. and Brosilow, C. (1987). Simulation of large scale dynamic systems - I . Modular integration methods. Compo Chern. Eng., 11 (3):241-253. Pantelides, C. (1988). Speedup - Recent advances in process simulation. Comp. Chem. Eng., 12 (7):745755. Perkins, J. and Sargent, R. (1982). Speedup : a computer program for steady-state and dynamic simulation and design of chemical processes. AIChE Symposium Series, 78 (214):000-000. Petzold, 1. (1983). Dassl : differontial algebraic system solver. Technical report, Sandia National Laboratories, Livermore, California. Petzold, L. (1988). An adaptativc moving grid method for one dimensional systems of partial differential equations and its numerical solution. III WOl'k.~hop on Adaptatiue Methods for Partial Differential Bquatunis, Oct. 1988, Rensselaer Polyruchnic Institute. Piela, P., Epperly, T., Westerbl'fg. h~., aud Wl'stprherg, A. (1991). ASCEND : an object oriented computer environment for modeling and analysis: tho modeling language. Compo Chem. Eng., 15 (1):53-72. Ponton, J. and Vasek, V. (1986). A two-level approach to chouiical plant. and process simulation. Compo Chem. Eng., 10 (3):277-286. Prokopakis, G. J. and Seider, W. D. (1981). Adaptative semi-implicit runge-kutta met.hods for the solution of stiff ordinary differential equations. Ind. Eng. Cheui. Fund., 20(3 ):255-266. Seader, J. D. (1989). The rare-based approach for modeling staged separations. Chern. Eng. Prog., October:41-49. Smith, G. and Morton, \V. (1988). Dynamic simulation using all equation-oriented approach to process flowsheeting. Convp . Chem. Eng., 12 (5):469-473. Westerberg, A. and Benjamin, D. (1985). Thoughts on a future equation-oriented flowsheeting system. Compo Chem. Eng., 9 (5):517-526. Wolbert, D., Joulia, X., Koehrot., I3 .. aud Pons, M. (1991). Analyse de sensibilite pour l'optimisation de precedes chimiques, In lUc(~nt.~ P'/'Ogrl~s en G(h!'i(~ d/~s P'I'odd(;'~, Contviegne 1991, pages 415-420.