Global modular Newton—Raphson technique for simulation of an interconnected plant applied to complex rectification columns

Global modular Newton—Raphson technique for simulation of an interconnected plant applied to complex rectification columns

Chid EngineeringScience,1976,Vol. 31, pp. 277-W. PergamonPress. Printed in Great Britain GLOBAL MODULAR NEWTON-RAPHSON TECHNIQUE FOR SIMULATION OF A...

688KB Sizes 6 Downloads 27 Views

Chid

EngineeringScience,1976,Vol. 31, pp. 277-W. PergamonPress. Printed in Great Britain

GLOBAL MODULAR NEWTON-RAPHSON TECHNIQUE FOR SIMULATION OF AN INTERCONNECTED PLANT APPLIED TO COMPLEX RECTIFICATION COLUMNS M. KUBICEK, V. HLAVACEKt and F. PROCHASKA Department of Chemical Engineering, Institute of Chemical Technology, 16628Prague 6, Czechoslovakia (Receiued 24 February 1975;accepted in revised form 15 October 1975)

Abstract-The module-oriented approach to steady state process simulation making use of the global Newton-Raphson method is described. The Newton-Raphson technique gives rise to a flexible and reliable computer program which can be used to solve highly nonlinear heat and material balances. The simulation system consists of a library of subroutines which are capable of automatic generating the residual vector and Jacobian matrix as well. The resulting set of linear equations governing the Newton correction exhibits a sparse structure (sometimes almost band) if the equations are grouped by modules. A mod&d Gaussian elimination was developed to solve a band matrix having some off-band submatrices. The concept of the generalized Newton-Raphson approach was used to construct a general computer code to solve rectification problems. The important point is the possibility of a straightforward solution of the design problems or problems of “controlled simulation” which via the flowsheeting method cannot be handled easily.

1.INTRODUCTION Digitalcomputers

are now being used to solve a number of problems concerned with the design and operation of a chemical plant. To treat a complex process the modular approach is usually adopted, i.e. the process equations are grouped into sets which correspond to a particular unit operation[4]. Each of these unit operations is programmed as an independent subroutine which is sometimes referred to as the “building block”. Any large process could be simulated by assembling these building blocks in an appropriate order. In solving feedback problems (e.g. recycle problems, counter-current problems, etc.) a master program (also an executive program) calls the unit operation programs in a sequence specified. Generally speaking, this method, when programmed, gives rise to a code which is not flexible enough to perform all types of calculations which the designer is faced with. Actually, for this type of calculation the distribution between two categories, dependent (unknown) and independent (specified) variables are fixed. In addition, the programmed dependent and independent variables need not be, of course, the very set which is required to answer questions for a particular design. Hence, to get a reliable tool for solving complex heat and material balances it is necessary to take advantage of the global, e.g. Newton-Raphson technique which is capable of providing versatility of the calculation, unlike the other methods (as, e.g. successive substitution, false transient method, etc.) which cannot yield easily the solution for certain specification of dependent and independent variables. The goal of this paper is to demonstrate the advantages of the Newton-Raphson technique to solve in a flexible way a large set of algebraic nonlinear equations describing the steady-state heat and mass balances in a complex plant. To illustrate this concept the Newton-Raphson

technique will be used to handle various problems arising in the design of rectification towers. A generalization to an arbitrary configuration of unit operations will be presented in a future paper. 2.Dl?PENDENTANDlNDEPFNDENl'VARIARLRS

A complex plant, which is represented by a flow-sheet, consists of a certain number of nodes (unit operations) which are interlinked by oriented streams. Each particular stream is described by a vector of variables containing, e.g. flow rate, composition, enthalpy and temperature of the mixture, pressure, etc. The jth stream variables will be denoted X+ Furthermore, let us denote X a vector including all stream vectors along with the vectors of parameters of all unit operations A: x = (X,‘, x2=, . . * , x,‘, P,‘, . . . , PK=y= (x,, . . . , XN)T Here J and K are the number of streams and nodes, respectively. Of course, for some types of design problems all parameters PK need not be specified. To perform a design study a certain number of variables are specified (independent variables), the other variables are considered as unknown (dependent variables). For instance, in calculation of the simplest conliguration of a distillation tower the feed’parameters, side stream rates and heat duties are independent variables while stage concentrations, temperatures and loadings are dependent variables. Basically, the dependent and independent variables can be interchanged, however, this operation must be done with great care because for some set of independent variables the set of equations need not exhibit a solution. 3.SYSTEM OF LINEARIZEDEQUATIONS

tTo whom all correspondence should be sent.

Let the streams jl, jz, . . . , j,,, be co-incident with the node k. The node can be described by the following 211

M. KUBEEK et al.

218

a

transformation relations: KIX,,, Jr-,, . . . Ix,,

&I= 0.

(2)

Equation (2) constitutes a set of nonlinear equations for k = 1,2,. . . , K. After grouping eqn (2) and denoting F = (F,‘, F*T,. . . , FK=y

(3)

the set of nonlinear equations becomes F(X) = 0.

(4)

With regard to the topology of a complex plant under question we may note that the structure matrix corresponding to eqn (4) is highly sparse, for instance, for a complex column a band matrix results with a low number of off-band elements. Let dim X = N and dim F = M, generally speaking, N-M variables can be chosen, i.e. we have N-M ihdependent and M dependent variables. The selection of the dependent variables is not completely arbitrary. For instance, the selection need not lead to an elimination of an equation as a result of a selection of such a set of independent variables which eliminates all variables in the particular equation. Evidently for this selection of variables an underdetermined system results. Of course, the selection of N-M independent variables may give rise also to linearly dependent equations. Finally, an inadequate selection of independent variables can lead to systems which do not yield a unique (or any) solution. We are not going to discuss in detail all possibilities of selection of variables, however, this topic would merit attention in future investigation. To solve the set of nonlinear equations the NewtonRaphson technique will be used: I’,(X’) AX“ = - F(Xk).

(5)

Here I is a Jacobian matrix. To improve the convergence properties the calculated correction AX” can be mod&d Xk+‘=Xk+~AXk_

(6)

The value of A,A E (0, l), is chosen to produce a decrease in the norm, e.g. in the sum of squaresPI: II~(Xk+‘)ll<

IWk II.

(7)

Here the vector X is a M-dimensional vector of dependent variables, the independent variable are assumed to be given and are invariable from iteration to iteration. The Jacobian matrix IF is formed from the columns of the matrix a& 1x I

i=l,...,M j=l,...,N

pertaining to the dependent variables. 4. MODULAR APPROACH TO TBE CONSTRUCTION OF A GLOBAL SIMULATIONPROGRAhi

The flowsheeting programs which calculate the heat and mass balances in a complex plant by a sequential approach are based on a modular construction. The sequential modular approach is very simple, information is transferred from unit operation to unit operation and hence the master program may be easily constructed. It seems reasonable to retain the modular approach also for

global simulation program because of its flexibility. This approach will be referred to as the modular global approach. Below the typical features of a modular global approach will be presented. The master program is supplied with: (1) A list and specification of all modules and streams. Based on this specification the system variables x1,.x2,.. . , x, are ordered and coded. (2) Information on the dependent and independent variables. (3) An initial set of all variables. The independent variables are fixed while the dependent variables represent a first guess for the Newton-Raphson method. (4) The topology of the interconnection between the modules. The master program calls sequentially all modules in a sequence which may be prescribed. The particular module accomodates the values of all inlet and exit stream parameters as well as values of module parameters. Of course, the values of both dependent and independent variables are transmitted. The specific modul transfers to the master program two informations: (a) the residual value of E. (The number of these values is given by the number of equations describing the given module); (b) The values of non-zero elements in the Jacobian matrix {aF,/ax,}. The values of non-zero elements are stored in particular rows, of course, the derivatives of 6 with respect to variables which are not coincident with the module under question are omitted. The submatrix of the first derivatives which is transferred from the module to the main program is Mk X Nk. Here Mk is the number of equations describing the module, Nk is the number of variables which are coincident with the module. The master program completes the zero elements in the rows of the matrix {~F,/c?x~}.Obviously, this matrix is highly sparse since Nk Q N. Let us note that the module creates the partial derivatives {a&/aXj} with respect to all variables which are coincident with the given module, i.e. also the derivatives with respect to independent variables are produced. After checking all modules the master program performs a reduction of the matrix {8F,/8Xj}to IF. This is accomplished by elimination of the relevant columns. The master program is capable of calling a subroutine for solving a set of linear algebraic equations (eqn 5), and a new approximation of X is produced. The calculation is continued as long as the pre-assigned tolerance IJAR”II< 6 is not reached. Let us note that the logical operations between a particular module and the main program are repeated in each iteration step, evidently, the iteration process is hierarchically superior to logical algorithms coded in the executive program. Below a detailed description of the modular global simulation program for complex rectification columns is presented. 5. MODULESFOR DISTILLATIONPROBLRMS The

general considerations presented above will be exemplified in this and the following section on a complex rectification column. A rectification column is essentially

Global modular Newton-Raphson technique for simulation of an interconnected plant

219

composed of four types of modules: reboiler, partial or Jacobian matrix is strongly improved. This finding is in total condenser and equilibrium tray. Each of these an agreement with our decision to group the equations and modules produce a part of the residual vector F, as well as variables according to modules. To make the situation clearer consider a simple column consisting of a partial the rows of the matrix {aE/&,}. Each module is tested whether the particular tray is not co-incident with an -condenser, reboiler and two stages where a threeexternal stream (feed stream, side stream or recycle). In component ideal mixture is distilled. The feed is sent to the positive case this stream is included in the module stage 2. As we have already shown the equations description. The equations describing a module are describing each particular stage are grouped in a sequence: M,, M2, M,, EB, SF, ER while the variables grouped as follows: are stored in the following way Xi.19X;,, X,l, Li, &, I&.The component material balance Mi i=1,2,...,I extra variables--QC,, Q&, F,, XF,, . . ., XF3-constitute EB enthalpy balance the last columns in the structure matrix, see Fig. 1. For this SF liquid mole fraction definition particular configuration we have 30 variables whilst only 24 ER vapor-liquid equilibrium ratio equations can be formulated. Evidently, if the feed The equations coded in the aforementioned modules are parameters are specified two additional parameters must reported in Table 1. be given. The most simple specification is the case where both heat duties QC, and QR,aregiven. We can notice that for this particular instance a block-tridiagonal matrix 6. ORDERINGOF TER VARIABLESAND FQUATIONS FOR DISTILLATIONPROBLEMS results. However, agreat number of other specifications is Recently Naphtali and Sandholm[ l] have shown that if possible. Consider, for instance, a specification of certain the equations and variables are grouped according to plate variables which are contained in the band. The columns the economy of storage of non-zero elements in the which are specified are replaced by the corresponding Table 1. Equations describing rectification modules Total condenser M,,i = Vj+,Kj+Jj+,.i

- V&i - L,X,,,

i=l,...,I

(1)

-%= A ~~,ixi.iH,i -5 (v, +L,W,.,h,.i + Qc, i=l (-1 %=$,&-1

(2)

(3)

+=I = i: K,.,& - 1 ,=I

(4)

For the partial condenserthe term VJ,,, is replaced by V&,,,X,.i Reboiler 4.i = L,$-,,i

- L,.iX,.i - VjKi,&

J=,= i: Lx,-,,hj-,,i I=,

.

-i

i-l

i=l,...,I

(5)

&X,.,/I,.,

- j$ ~Kix,iH,i - QRi (6)

se=&,,-, i=, ERi = 2

,=I

K,+& - 1.

Equilibrium tray NJ= L,x,-8,

EBi = i

i-1

-65, + V&)X,.i

Li-1&-~.&-l.i -i

,=I

(L&i t vjK,JZ~,~)X~.i

+ i: V~+IK,+I.,H~+~‘C,+,., 1-1 SF,=&,*4 i-1 ER, = i: K,.,X,., - 1.

i-1

t y+,K,+,.iX,+,,i

(8)

i=l ,...,I

(7)

M. KL&EK et al. Off-band Band variables

xxxn

xxxn xxxn xxxn

xxx xxx xxxx xxxxxx xxx xxx

xx xx xx xx

X

X

X

X X

xxx X

xxx xxx xxxx xxxxxx X

X

x

xxx

xxx

X

Fig. 1. Basicconfigurationof the structurematrixcorrespondingto Jacobianmatrix(singlecolumn).

columns from the right-hand column submatrix. For example, if XI,*and L, are specified then after rearranging following matrix results (see Fig. 2). We may note that the sequence of variables is changed. The interchange of columns is undertaken in order to minimize the number of off-band elements. For a structure involving a sequence of complex interlinked rectification columns the appropriate grouping of variables may be, e.g. X ,,,, . . .)X,.1, Ll, VI, l-1, X2.1,* . . 9X 2.11

4 x & MS El3

x x

x X

xxx

xxx xxx xxxx xxxxxx xxx xxx

X

x

x

...,QRB, REC,, . QR ,,

,...

e

.

, REG.

X

if

xx

X

X x

xxx xxxx

xxx xxx

Ml

X

X

xxx

X

X

x

X

xxxx

xxx xxx .xxxx xxxxxx xxx xxx X

xx

MS

Xt

xx

X

x

Ma

X

xxxxxx

ER

xx xx

x

xxx X

xx xxxx

xx xx xx

xxx

x

x

,FE,

The reader may notice that this grouping is not optimal with regard to the occurence of off-band elements in the Jacobian matrix, for huge problems the off-band elements should be grouped according to columns.

X

X

X

Mt

ER

S&F,

XF ,,,, . . . , XF,,r, . . . , X&I, . . . , XFm,

x

X

Mn MS

z

QC,,SS I,...,

x

X x

QC ,,...,

xx xx xx xx

X

x

& MS ES SF ER

X

xxx

SF x ERX MI

xx xx xx xx

Lz, Vz, Tz, . . . , XJJ, . . . , -%I, LJ, V, TJ,

x

X

Fig. 2. Configurationof the modified structure matrix corresponding to Jacobian matrix (singlecolumn).

Global modular Newton-Raphson technique for simulation d an &?r&CWfiected pfant 7.

METHODOF !3OLUTION OF LINEAR ALMOSTBANDSYSFEMS

The algorithm is based on the technique of modified matrices. Let us have a set of linear equations

\

fx = b

R = R,RzT here R,,R2 are n x m matrices. The matrix R, is composed of non-zero columns j,, h,. . . , j,,, of f - T. The matrix R2 is formed from the unit vectors ei,, eh, . . . , % The algorithm performs the following steps: (1) The matrix V(n x m) and the vector y satisfy TV= R,

(10)

Ty= b.

(11)

The modiied Thomas-Gauss algorithm is used to split T=LU, V and y are calculated by back solving m t 1 times. In principle, this algorithm is a standard LU factorization of a multidiagonal matrix[7]. (2) Form a matrix A(m x m) (12)

and a vector w = Rz=y. (3) Solve AZ = w

for z.

(13)

(4) Calculate x = y - vz.

(14)

The method described is very powerful for m Q IZ,it can be used also if m
complex tower and interconnected columns as well. A detailed discussion of various engineering aspects of calculation of dilhcult separation problems (superfractionators, highly non-ideal mixtures, non-standard specifications, etc.) will be published elsewhere191.

(9)

where f = T + R, here T is a multidiagonal (band) matrix n x n and R is a matrix of low rank. Let us write

A =ZtR,=V

281

Example 1 Example 1 from Hanson’s book[3] has been chosen, i.e. a six component mixture is rectified in a column having 8 trays with a reboiler and a partial condenser. The feed is boiling liquid and is sent to stage 5. A side-stream is withdrawn from stage 3. Details on the vapor-liquid equilibrium constants are presented in the Hanson’s book. In this particular instance the following variables are specified (independent variables): V,, L1, SS, F, XF,, , . . , XF,. Of course, for this specification the relaxation method (i.e. the false transient method) can be easily used1121. The problem requires to solve 90 nonlinear algebraic equations. The structure matrix which is commensurate with the Jacobian matrix is similar to that drawn in Fig. 1. Evidently, here a six component mixture and a column having eight trays is considered. In addition, in contrast to the column specification displayed in Fig. 1, here a side-stream is given, i.e. the variables in the off-band column matrix are grouped as below: QC, SS, QR, F, XF,, . . . , XF,. With regard to the specification only the elements occuring in the QC and QR columns are unknown. These columns, of course, take the position of the V, and L, columns since the variables V, and L, are specified. As the reader may infer from Fig. 1 the unknown variable QR after the column interchange falls outside the non-zero band, similarly, as we have shown in Fig. 2. As the initial concentration profile the feed composition have been used at each stage while as an initial profile of temperature the linear profile has been guessed between 117 and 235°F. The initial flow rates have been calculated from specification. The solution has been reached in 3 iterations (CPU 50 set on IBM 370/145).After 3 iterations the value of the sum-square residuum decreased from 0.14 X 10’ to 0.6~ 10m5.On the other hand, to compute the solution with the same accuracy, the relaxation method requires 17 iterations (CPU 35 set on IBM 370/145). Example 2 This problem possesses a different specification. Here the pertinent indepent variables are: L1, X1,4, SS, F, XF,,..., XF,. The initial guesses are the same as in Example 1. Evidently, this problem cannot be solved in a straightforward way by the relaxation procedure. Gf course, since this specification is of the “controlled simulation” type the conventional procedures has to be combined with some interpolation technique which gives rise to lengthy calculations. The global approach, however, possesses no dithculties. To reach the solution within the tolerance (of composition) 0.001 seven iterations are sufficient. The CPU (IBM 370/145) has been 112sec.

282

M.

KUB@EK et al.

Example 3 In two interconnected distillation columns a six component mixture is fractionated. The mixture to be rectified is the model mixture presented in the Hanson’s book. Column 1 has a reboiler, 8 plates and a total condenser, whereas the column 2 is provided with a reboiler, 6 plates and a total condenser. The feed is fed as boiling liquid to the plate 5 in the column 1. The feed rate is 1 kmol/hr and the composition is XFI = Xl% = 0.1, XF, = XF, = XF, = XF, = 0.2. The bottom product is withdrawn and sent to the plate 3 of column 2. Liquid from the reboiler of column 2 is returned to stage 1 in column 1 and is set at 0.1 kmol/hr. The amount of the top products in columns 1 and 2 are 0.7 kmol/hr and 0.2 kmol/hr, respectively. For the first trial it is assumed that the concentrations at all stages are equal to the feed composition, the linear temperature profiles are guessed between terminal values 70 and 230°F as well as 215 and 245°F for the first and second columns, respectively. The initial flow rates have been calculated from ‘specification. A sketch of the configuration is depicted in Fig. 3. The underlaying specification is: L,, VI, L,,, VI,, F, XF,, . . . , XFS, REC,, REC,. For the purposes to acquaint the reader with the structure of the pertinent Jacobian matrix Fig. 4 is presented. Figure 4 reveals that the band-width is occupied by 27 elements, the occurence of the off-diagonal elements is obvious from Fig. 4. We have chosen this particular case because this configuration was calculated by the relaxation procedure earlier[lO] and it is a possibility of comparison. While the relaxation method yields convergent solution in approx. 25 iterations, the Qc

Newton-Raphson technique is capable of getting the solution in 7 iterations. Of course, for specifications which, for instance, 6x the concentration of one component in the exit streams the conventional methods cannot be economically used. A detailed discussion of this problem is presented elsewhere[9]. 9. CONCLUSIONSAND DISCUSSION

Whereas the computer-aided-design literature abounds with the modular flowsheeting approach to the numerical solution of simulation problems, so far relatively little guidance is available in the literature that may be useful in adopting the Newtondaphson approach to simulation of an interconnected plant. Most of the information has been acquired for the calculation of single complex rectification towers[l, 21. Current acceptance of the global modular Newton-Raphson approach depends to a large extent on following items: (1) development of powerful sparse matrix procedures. (2) possibility of non-numerical mathematical operations and automatic generating of analytical expressions for the elements of the Jacobian matrix via algebra packages[ll]. (3) use of a computer with a large core capacity. The principal considerations in choosing the flowsheeting approach are simplicity, modularity, flexibility and convergence. However, this method has the following deficiences: (1) For physical quantities which are dependent on all stream parameters (e.g. nonideal vapor-liquid equilibrium) there is no evidence of convergence. (2) It is not feasible to use this approach to “controlled simulation” problems where more than two stream parameters must be adjusted. (3) Sometimes the convergence is very slow. On the other hand, the Newton-Raphson approach has the following virtues: (1) It is computationally advantageous since each admissible specification can be calculated. (2) The rate of convergence is high. However, much remain to be done. For instance, there is an urgent need to extend the Newton-Raphson approach to such equations which cannot be easily differentiated, e.g. to highly nonideal rectification or extraction problems or problems which embrace in the description diierential equations, the theory allowing to classify the specification ((a) underdetermined (b) overdetermined (c) correctly determined) of a particular problem is missing as yet, etc. There are still obstacles to be overcome before this approach can be virtually a commonplace item in the modem computer-aided-design. Nonetheless, the authors believe that it is undoubtedly possible and that it merits careful study in future. NOTATION

Fig.3.A sketchof twointerlinkedrectificationcolumns.

molar flow rate of feed i molar enthalpy of component i in liquid at stage j molar enthalpy of component i in vapor at stage j equilibrium constant of component i at stage j’ molar flow rate of liquid phase at stage j number of equations

Global modular Newton-Raphson technique for simulation of an interconnected plant Condenser Stage1 2 3 4 5 6 7 6 Reboile? stage 1 2 3 4 f Reboiler

Submatrix As

Submatrix AP 5

-. 7 s

diTsc

.

X X X X X

xxx xxx xxx xxx xxx xxxx

xxxxxxxxx xxxxxx xxxxxx

xx xx xx xx xx xx xx

X X X X

X X

xxxxxx X

gr

i;

. . . . *‘JSC

X

X

X

X X

X

X

X

x

x xx

xxxxxxx

I

x

L

Submatrix A4

Submatrix As

SubmatrixA.

SJSC‘

.

G.. X

X

X X

X X

X

X

X

X

X

xxxxxx

X

xxxxxx

SubmatrixA,

Fig.4. Conligurationof themodiiiedstructurematrixcorrespondingto Jacobianmatrix(twointerconnectedcolumns).

N

number of variables vector of parameters of unit operation i heat duty in condenser i QRi heat duty in reboiler i RECi molar flow rate of recycle stream i Ss;: molar flow rate of side-stream i vj molar flow rate of vapor phase at stage j Xi ith component of vector X, see eqn (1) x vector including stream vectors along with vectors of parameters of all unit operations vector of stream i molar fraction of component i at stage j in liquid phase

Q2

molar fraction of component i in feed j tolerance A coefficient in eqn (6) rp Jacobian matrix

X&

l

REFERENCES

[l] Naphtali L. M. and Sandhohn D. P., A.I.C.H.E. J. 1971 17 148. [2] Goldstein R. P. and Stanfield R. B., Ind Engng Chem. Proc. Des. Dev. 19709 78. [31 Hanson D. N., Somerville G. F. and DuKin J. H., Computation of Multistage Separation Processes. Reinhold, New York 1962. [4] Crowe C. M. et al., ChemicalPlant Simulation. Prentice Hall Englewood Cliffs, N.J. 1971.

284

M. KUBEEKet al.

[S] Thomas L. H., Dept. of Watson Sci. ComputingLab., New York 1949. [6] Forsythe G. E. and Moler C. B., Computer Solution ofLinear Algebraic Systems. Prentice Hall, Englewood Cliffs, NJ. [7] kz&kM Algorithm470 Communs ACM 197316,12,760. [8] KubfEek MI, Communs A&f, to be published. [9] Prochlska F., Eckert E., Hlav%ek V. and KubfEekM., Ind.

Engng Chem. Proc. Des. Dev., in preparation. [10] Jelfnek J., HlavQek V. and Kfivsky Z., Chem. Engng Sci. 197328 1833. [ll] Xenakis J., IBM Contributed Program Library. IBM type III. No 360 D-03.3.004 (IBM Boston Programming Center, Cambridge, Massachusetts). [12] Jelfnek J., HlavYek V. and Kubfc’ekM., Chem. Engng. Sci. 197323 1825.