Computers & Chemical Engineering, Vol. 3, pp. 535-542, Printed in Great Britain. All rights reserved.

1979

0098&1354/79/040535-08$02.0/O Copyright @ 1981 Pergamon Press Ltd.

Paper 10B.l NAPHTALI-SANDHOLM DISTILLATION CALCULATIONS FOR NGL MIXTURES NEAR THE CRITICAL REGION LARS J. CHRISTIANSEN Haldor TopsBe A/S, Nymcallevej 55, DK-2800 Lyngby,

Denmark

and MICHAEL L. MICHELSENand AAGE FREDENSLUND* Instituttet

for Kemiteknik,

Danmarks

Tekniske

Hsjskole.

DK-2800

Lyngby,

Denmark

(Received 2 July 1979)

Abstract-Multicomponent distillation calculations are carried out using a modification of the Naphtali-Sandholm column calculation procedure coupled with the Soave-Redlich-Kwong equation of state for generating K-values and enthalpies. The modified Naphtali-Sandholm procedure is applicable when the individual stage efficiencies are unity and the K-values are not strongly dependent on the vapour-phase composition. For a column with M components and N stages, N(M + 2) relationships are set up. The N(M + 2) equations are solved using Newton-Raphson iteration with imposed maximum allowable changes of the independent variables between iterations. Analytical derivatives of the K-values and enthalpies with respect to composition and temperature are employed in the Jacobian. The method has been used for the process development and design of distillation columns associated with the separation of natural gas liquids (NGL). Fast and safe convergence is obtained even near the critical region. Scope-Existing algorithms for solving multistage, multicomponent separation problems can be divided into three main groups, each with its range of applications. The algorithms differ in the methods used for converging the material and energy balances and the equilibrium relationships. The socalled O-method of Holland et al. [ 11 uses sequential convergence. From an assumed set of stage temperatures and interstage flows, equilibrium ratios on all stages for all components are calculated. Next, the component material balances are closed one by one by the tridiagonal Thomas algorithm, and the improved compositions thus obtained are used to correct the stage temperatures. Finally, the total flows are corrected using the enthalpy balances. The whole process is repeated until convergence is obtained. Convergence is fast for conventional distillation of fairly narrow boiling, nearly ideal mixtures. An excessively large number of iterations may be required for nonideal mixtures, wide boiling mixtures, and for absorber and stripper calculations. The simultaneous convergence method of Tierney and Yasonik [2] is applicable to all types of separation processes and has very good convergence properties for nearly ideal mixtures. The stage temperatures and the total liquid flows are taken as independent variables giving a total of 2N independent variables. The 2N check functions used are the sum of vapour phase mole fractions and the discrepancy in the enthalpy balance for each stage. The equations are solved using Newton-Raphson iteration. The method is particularly well suited for separation units with fairly few stages and (nearly) ideal mixtures, for which convergence is quadratic. The most general approach by far is that of Naphtali and Sandholm [3]. Equations are grouped by stages rather than by components and solved by the Newton-Raphson method. The number of independent variables is very large, for a system containing M components it is N. (2M + 1). The resulting matrix of partial derivatives is block-tridiagonal, with blocks of size (2M + 1). (2M + 1). The method is applicable to any kind of countercurrent, stagewise separation, and the order of the convergence is not affected by nonideality. Contrary to the other methods, evaluation of partial derivatives of all thermodynamic properties with respect to the independent variables are required, a handicap which in our opinion is more apparent than real. The method was originally intended for highly nonideal mixtures. However, computation time and storage requirements depend much more strongly on the number of components M than on the number of stages N, and hence the method may be superior to the other methods even for relatively ideal mixtures when N is large and M small. Frequently substantial simplifications of the Naphtali-Sandholm method are possible. A simplified version assuming constant molal overflow and ideal stages is described by Fredenslund et al. [4]. Gallun et al. [5] describe an interesting variant, the “almost band algorithm”, applicable for columns without side streams. Total material balances around one *Author to whom correspondence CACE3: 1-4 - II

should be addressed. 535

536

LARS J. CHRISTIANSEN

et al.

end of the column are used to eliminate the vapour phase component flows. This reduces the block size to (M + 2). (A4 + 2). A block size reduction (M + 2). (A4 + 2), equal in size to that of Gallun et al. but applicable to columns with arbitrary side streams, is described in this work. Ideal stages are assumed. The convergence properties are only adversely affected by the simplifications in the rather unusual case when the equilibrium ratios are strongly dependent on the vapour phase composition. Conclusions and Significance-A modification of the Naphtali-Sandholm method for stagewise, countercurrent column calculations has been proposed. The modification is intended to be used when the K-values are relatively independent of the vapour phase composition and the stages are ideal. The computer time requirements are about a of those for the full NaphtaliL3andholm method. The Soave-RedlichhKwong equation of state is written in a form especially suitable for computer applications. Analytical derivatives of K-values and enthalpy-deviations with respect to temperature and composition are employed. This has been found to greatly facilitate the column calculations. Similar experience is reported elsewhere for the UNIQUAC equation [l 11, and the method may be extended to other equations (Peng-Robinson, Wilson, etc.). Finally, practical application of the algorithms to the design of natural gas liquid plants is illustrated. COLlJMN

MATERIAL-AND

be distillation other

column,

ENERGY

N

BALANCES

DISTILLATE

CCNDENSER

in Fig. 1 or any

a be introduced

from

at be withdrawn or cooling

i the independent M liquid phase component 1, the stage temperature T, and the total vapour phase flow rate y. The nomenclature for a typical stage is shown in Fig. 2. The following balance equations are set up : M component material balances

FEED

SV nc1

nt1 L

FEED

n

Fn'fn,i

G-1 n-1

:

ZJ lij(l

+~)+~ij(l+~)

-li+,,j

-

ui-,,j

-Jj

=

The equilibrium relationship vi _ , .j (see Eq. 5).

is substituted

One bubble point relationship

:

c

Ki,

li,

-

c

0.

(1)

for uij and

BOTTOMS Ii,

=

0.

(2)

k

k

Fig. 1. Distillation

Subscript

i:

Subscript

j: component

flow

from

stage

column

i, i =

1, 2.

j, j = 1, 2,.

V = total vapor flow v = component

I

L = total liquid

1 Stage

I = component

I i

F = total

k-f,

1.

.sl’

I

L

h., hn

vapor

flow

flow liquid

flow

feed

SL = liquid

side stream

S” = vapor

side stream

h = liquid enthalpy H = vapor enthalpy

I

Vi-, vi- 1 ,j H,-,

Fig. 2. Nomenclature

for an arbitrary

stage in a distillation

column.

M

configuration.

N

PRODUCT

Naphtali-Sandholm

distillation

calculations

One enthalpy halunce: (Li+SLi)hi+(~+Sy)Hi - l’_,Hi-l

-Li+lki+*

- Qi = 0.

(3)

Li is the total liquid flow rate, and SL, and SF denote liquid and vapour side streams from stage i. Fi andfij are total and component feed rates to the stage, and Qi represents the enthalpy added to the stage (with the feed, by interstage cooling, or otherwise). The dependent variables Li and vij are calculated from

for NGL mixtures near the critical region

For a column with 25 stages and 9 components, convergence to relative changes in the independent variables of less than 10d6 normally requires 336 iterations and about 4 sec. on the IBM 370/165. The program requires 70 K bytes storage. Naturally, neglecting the partial derivatives of the K-values with respect to the vapour phase composition does affect the rate of convergence. However, the reduction in dimensionality provides ample compensation. The present modification is typically four times faster than the full NaphtaliiSandholm procedure. OF THE SOAVE-REDLICH-KWONG

APPLICATION

Li =

oij =

qyij

=

EQUATION

1 lij

v~ 1

Kijl,,

(5)

Kiklik

k

where Kij is the equilibrium ratio y/x for component j on the ith stage. Calculations are normally performed with specification of all feeds, all heat inputs, and all product including the distillate and bottoms. streams, However, the energy balances around the bottom and top stages are replaced by the following relationships : C /,,j = L, i

(bottom

product)

(6)

T I,, = B(F,V + SLV) where R is the reflux ratio. The N(M + 2) equations corresponding to Eqs. (l)(3) are solved by NewtonRaphson iteration. Evaluation of the Jacobian matrix requires computation of partial derivatives of these relationships with respect to the independent variables. This means, among others, that a method for evaluating derivatives of equilibrium ratios and enthalpies with respect to composition and temperature must be available (see Appendix B). Solving the N(M + 2) equations by NewtonRaphson iteration is not as formidable a task as it at first might appear to be. We notice that the equations on stage i are only affected by the variables on stages i - 1, i and i + 1. Hence the Jacobian matrix of partial derivatives becomes block tridiagonal with block size (M + 2). (M + 2). The partial derivatives are here found by differentiation of Eqs. (1 t(7) employing analytical derivatives of K,,, Hi and hi with respect to the independent variables. Partial derivatives of Kij with respect to yik are ignored. In the computation of K-values, y-values from the previous iteration after a single correction step are used. Initial estimates of the total interstage flows are found from the total material balances, assuming constant molal overflow, and a linear temperature profile is imposed. K-values based on the composition of the combined feeds are used for solving the component balances in the first iteration step using the Thomas method. The initial temperature and composition profile obtained in this manner can normally be applied directly to the modified Naphtalii Sandholm procedure.

537

OF STATE

Equilibrium ratios Ki, and enthalpies Hi and hi are needed in the solution of Eqs. (1 t(7). We choose to use the Soave-Redlich-Kwong equation of state [6] for evaluation of the equilibrium ratios and the enthalpydepartures from ideal gas-values. Two-constant equations of state such as the Soave modification of the Redlich-Kwong equation can not be highly accurate in the whole range of temperature, pressure, and composition. But they have the advantage that they require relatively little computer time, and the Soave-Redlich-Kwong equation seldom gives rise to large errors when it is applied to mixtures of lower hydrocarbons, nitrogen, carbon monoxide, carbon dioxide, and hydrogen sulphide. The same model is used for the vapour and the liquid phases. Thus-unlike for example with the ChaoSeader method-continuity at the critical point is retained. The SoaveeRedlich-Kwong equation is stated in a form particularly useful for computer applications below. For a mixture, the equation of state is : PV

V ~V-b

‘=R,T=

(ad R,T(V+h)

(8)

where - kkj)

(a~) = c 1 xjx,aja,zja,(l j k h = 1 xjhj; j

hj = 0.08664 RgT,j/P,j

aj = CO.42747 R: ‘i$/P,J’.’ aj=

l+~~(l-T,qi.~);

Tj=

T/T,,

and m, = 0.48 + 1.574~~ - 0.176wf. The necessary critical constants and binary interaction parameters kij are given by Reid et al. [7]. For hydrocarbon-hydrocarbon interactions, kij is zero. The double sum (aa) is time-consuming to evaluate on the computer. It is reduced to a single sum as follows : (aa) = 1 xi*zT where .xi* = xjajaj and zf = C x:(1 - kjk). L

(9)

538

LARS J. CHRISTIANSEN

o~C

~

et al.

--

~ u

Z to) t x l t N

t..)

~" (' 4 (' 4

U') n~

O

N

N

u.I 121

~

O

O t~

z

<

~

°

~

o

Naphtali-Sandholm

distillation

calculations

The figacity coefficient is : In $j = i

(2 - 1) - ln(z - B) -

i[v-F]ln(l

+s).

(10)

The enthalpy departure (vapour and liquid) is : AH

_--_~~_~--_ R,T

1 bR,T

xcx;zj* l+mj. c(,

( >

i

(11)

J

In the above two equations A = (aa)P/Ri TZ The K-values lated from :

and

and enthalpies

B = bP/R,T. on stage i are calcu-

Kij = f$$$;

(12)

Hi = Ho - AHi

(13)

hi = Ho - Ahi

(14)

where Ho is the ideal gas enthalpy. Expressions for the derivatives of the enthalpy departures and K-values with respect to temperature and liquid phase composition are given in Appendix A. We use analytical differentiation rather than numerical for several reasons. The analytical derivatives are accurate, their evaluation is much less time consuming on the computer, and convergence near the critical point of the mixture is not adversely affected by numerical manipulations. Evaluation of all fugacity coefficients and enthalpies for a nine component mixture requires 0.7 msec. on the IBM 370/165. Evaluation of all partial derivatives requires additionally 1.5 msec. compared to 7 msec. for numerical differentiation. For the program used in this work it is estimated that only about 25% of the computer time in the distillation calculations is spent on the evaluation of thermodynamic properties. The Soave-Redlich-Kwong equation of state has the property that the critical compressibility factor for a pure component is 3. This is not necessarily the case for mixtures (nor for many pure components), but since distillation columns in practice are not operated very close to the critical point, where separation is adversely affected, we may safely assume, that a liquid phase exists if the corresponding calculated compressibility factor is less than 3 and that a vapour phase exists ifit is larger than f If the calculations during the column iterations have resulted in a phase which does not meet these requirements, the phase is hypothetical and the temperature is changed to a level where both phases exist. The above procedure ensures that reliable thermodynamic properties for real (non-hypothetical) phases are found efficiently and that the column calculations will converge even rather close to the critical point. For example, the depropanizer of Fig. 3 has been converged about 5 atm from the critical pressure of the distillate. NGL SEPARATIONS

A successful, practical application of the algorithms described above is illustrated in this section. The

for NGL mixtures

near the critical region

539

application concerns the design of a Natural Gas Liquid (NGL) plant, and Fig. 3 shows the plant configuration chosen for illustration. The feed is treated natural gas containing saturated hydrocarbons ranging from methane to heptane. A gas containing mostly methane and some ethane is separated overhead in the demethanizer; the demethanizer bottoms product is separated so that ethane and propane goes overhead. Consequently, the demethanizer yields the ethane and propane products; the debutanizer yields the butane product in the distillate, and the bottoms product is the heavy ‘gasoline fraction’. To avoid excessive refrigeration costs, the pressures are everywhere, except for the debutanizer, relatively high. This means that most of the separations are carried out near the critical pressure. The design problem is now for each column to find the number of stages, reflux ratio, feed stage location, and total product split so that the specified products may be obtained most economically. Finally, reasonable column pressures, governed by the mixture thermodynamics and economics, must be established. The revised Naphtali-Sandholm procedure has been included in a program SEPTOWER [S, 93, which among other methods can utilize the Soave-RedlichKwong thermodynamics. The program calculates the performance of a single distillation column including the reboiler and the condenser. The individual distillation columns are coupled together by use of a computer operating system called GIPS, General Integrated Programming System [lo]. This system permits the user to build up a calculation from an arbitrary combination of existing programs, and it uses a common data base for all programs which makes it easy to transfer data from one program to another during a calculation. The system also contains generalized input-output facilities. The GIPS system is used to combine the columns shown in Fig. 2. The product from one column is by simple instructions used as input to the next column. The state of any feed or product stream may be changed from one temperature, pressure, and enthalpy to another. The latter requires the use of a flash program, but this is easily included in the calculations due to the flexibility of GIPS. The column design variables: number of stages, reflux ratio, feed stage, and total product split are assigned by trial and error to fulfil the product requirements. NOMENCLATURE

SRK parameter mixture SRK parameter (aa) SRK parameter b SRK parameter I3 mixture SRK parameter component feed f F total feed on stage H vapour molar enthalpy h liquid molar enthalpy AH enthalpy departure from ideal gas Ho ideal gas enthalpy K K-value, K = JJ/X I liquid component molar flow from a stage L total liquid molar flow from a stage SRK parameter .G number of components N number of stages

A”

LARS J. CHRISTIANSENetai

540

where

pRSS”*e

i R, SL sv T T, t: V x Y z

interstage heating or cooling reflux ratio gas constant liquid side stream vapour side stream temperature reduced temperature, T, = T/T, component vapour molar flow total vapour molar flow molar volume liquid mole fraction vapour mole fraction compressibility factor

zxll

zx12 = 1 x,jzrj;

= ~Xljz,,; J

This allows the enthalpy AH --Z=sRT

zx22 = c XZjZZj.

I

J departure

l -&[ln(l

to be written as

+!$][zxll

:

+JT.z.x12).

(44 Derivatives of SRK Derivatives of 4i and AH with respect to temperature and composition are needed in column and flash calculations. The following derivatives are needed later :

Greek letters a: SRK-parameter w acentric factor 4 fugacity coefficient

c’b __b;; (7.X,

Subscripts c critical i stage number j, h component number r reduced

(54

!!!=!?&g; ?X;

t

(64

c;zc?A CiZ -=zz++‘B:?( SXi '_I

8; irB '

(74

L

Since Eq. (8) of the main article may be written as REFERENCES

5.

6. 7.

8. 9. 10.

Il.

C. D. Holland, Fundamentals and Modeling of Separation Processes. Prentice Hall, New Jersey (1975). J. W. Tierney & J. L. Yasonik, AIChE J. 16,897 (1969) L. M. Naphtali & D. P. Sandholm, AIChE J. 17, 148 (1971) Aa. Fredenslund, J. Gmehling, M. L. Michelsen, P. Rasmussen & J. M. Prausnitz, IEC Process Design and Development 16,450 (1977). S. E. Gallun, F. E. Hess, G. W. Entzen & C. D. Holland, 5th Symp. Comput. Chem. Engng, High Tatras, Czechoslovakia, Octover 1977. G. Soave, Chem. Engng Sci. 27, 1197 (1972). R. C. Reid, J. M. Prausnitz & T. K. Sherwook, The Properties qf Gases und Liquids, 3rd Edn. McGraw-Hill, New York (1977). J. K&r, Chemical Engineering Computer Program. Haldor Topsoe, Ravnholm, Denmark (1975). Topsee Topics, August (1978). J. Kjar & P. Sorensen, Introduction to the General Integrated Programming System. GIPS (2nd Edn), Haldor Topsae, Vedbsek. Denmark (1969). T. Magnussen, M. L. Michelsen & Aa. Fredenslund. Distillation Symp. ‘79, The Institution of Chemical Engineers (London) Symp. Series 56 4.2/l (1979).

it follows that B-z c’z -=__’ 8A N

(;7z -=

A+z(1+2B)

’ 8B

where N=z(3z-2)+A-B(l+B) c:B

B

f?T

T

__(7;;) - 2x22 + zx 12/J? dA r=A

&-l [

1 ;

Q4=

-zx114%l2

Temperature derivative of thefigacity coeficient From Eq. (10) of the main article it follows that

(154

Equation (8) gives the SRK equation of state for a mixture, and the fugacity coefficients and enthalpy departures are given by Eqs. (10) and (1 l), respectively. The term (a~() given by Eq. (9) may also be expressed as follows : Let for pure component i = ali + u,~JF

Pa)

N

APPENDIX A THE SRK EQUATION OF STATE

(UC+ = a,(1 + mi)-(a,m,T,,0.5)fi

(84

z3 - zZ + (A - B - B’)z - Al3 = 0

Only the term written as

DA needs

further

evaluation.

DA =QIQ;+Q'IQz

DA may

be

(16a)

where

(la) (l7a)

If XT =

ctliXi

+

(U2J,)~ JT

= x,i

+

xz;fi

(184

and $=CxIj(l-kij)+ j (24 then it may be shown that (aa) = zxll

+ 2fi.zx12

+ T.zx22

(3a)

NaphtaliiSandholm

distillation

calculations

Temperature derivative of the enthalpy departure This derivative is evaluated from Eq. (4a) :

-

g+R ( T$’

I

-

>

-;QI(ax12,,i7;)

k Q’,(zxll

+ z.xl2~fi).

Composition derivatives ofthefugacity d In 4i 3Xj

1

bi

dz = [p-%-1) b dxj b

(22a)

coeficienr

-Qs-Q&rQkQ, 6%

for NGL mixtures

The Jacobi matrix. (ZF/Ex), is here very large, but its evaluation is greatly facilitated by the fact that the conditions on stage i only are influenced directly by the conditions on stages i+ 1 and i- I. As a result, the Jacobian becomes block-tridiagonal in structure, which permits rapid solution by block elimination. The diagonal elements of the Jacobian, B, contain the derivatives for stage i with respect to the variables on stage i. The elements below the diagonal, A, contain the derivatives for stage i with respect to the variables on stage i - 1. The elements above the diagonal, C, contain the derivatives for stage i with respect to the variables on stage i + 1. The structure ofthe Jacobian and the system ofequations is then :

B,C, AAC, A,B,CS ‘.

where

P-W

A,_,B,_,C,_

(254 (264

(274 In K, = In f#$ - In 4:. In the SRK distillation ation is made :

program,

approxim-

5 In 4,”

~3In K, axi

(2ga)

the following

-

(29a)

dx,

Composition derivatives of the enthalpy departure

-

F =RT 2 - QsQ, - QLQ, J

(304

J

where

(314 Qk = -

bI 2 Qa + P(a-a-z h J I IJ +

h

Partial molar quantities are calculated

SOLUTION

ccljzT),

according

541

near the critical region

I

A,

D, D2 4

1* I

D,

1

Bw

D,

ADi are the increments in the independent variables, x’“~“” - x”“~‘, for stage i, and Di are the values of the function Fi (ideally zero). The individual elements of the Jacobian are explained below. Theequations are solved by Thomaselimination. Row N is multiplied by EN ’ whereby the diagonal matrix B, is reduced to the identity matrix. The resulting row is multiplied by - C,V_ ,, and this result is subsequently used to eliminate C,_, from row N - 1. The procedure continues to row 1, from which C, is eliminated yielding directly a value for AD,. The remaining increments AD, are formed by back substitution from row I to row N. The increments AD, are added to the old values of the independent variables. The iteration continues until Zi AD; is smaller than a given tolerance. Derivation ofthe Jacobian (all elements not shown are zero) For each stage, the first M equations are the component mass balances; then follows the bubble point relationship and the energy balance. The first M independent variables are the liquid component flow rates; then follows the temperature and the total vapor flow rate. In the following, Ml = M + 1, M2 = M + 2, i is the stage number, j the equation number, and k is the index for the independent variables. Material balance (Eq. 1) i = I, 2.. N

(324

to Eq. (19b).

B,, = 2

,k

(v, + SK) + 6,( 1 + SL,/L,) - l,,SL,,‘Lf

I’?~,is the Kronecker

APPENDIX B. OF THE COLUMN DESIGN EQUATIONS

(4’4 (%I

(lb)

where F is the vector of functions (I)(3). In Newton-Raphson iteration, a new set of values of functions Ft”‘“’ are generated from a previous estimate in the following fashion

(new,_ xi”‘d)) = 0, (x Yi=XIYld>

(W

delta.

Solving Eqs. (l)(3) of the main article means finding the set of independent variables x(lij, 7; and I’,) which satisfies the relationships F(x) = 0

W

(2b)

This equation is used to estimate x”‘~~‘. When xrneWJ- x”“~’ is sufficiently small, the correct set of values of x has been found, and the iteration stops. In the computer program, x’“~‘“)- x”ld) is, between iterations, limited in the following way: The temperature is allowed to change max. lO”C, and the liquid component flow ratio are allowed to change max. 20 90.

(W c,,, = -6,. Bubble point relationship

Vb)

(Eq. 2) i = I, 2

Bi.ncI.k = c l,j 2 I

N

+ K,, - 1;

B ,,u,,u,=I:I,,s!.

KW

I Energy balance

(Eq. 3) i = 2.3,.

B I.MZ.I= f$

,I

N - 1

^ (v, + SF) + $(Li

d

+ SL,) + h,

W-9

542

LARS J. CHRISTIANSENet al.

sYij 7&f = Yij

(lob)

AI.MZ,MZ = -H,-, ; &,Z,MZ = Hi A I.h42,k = -vi_, A uuZ.MI=

_v

(1lb)

The derivatives

UW

of the vapor enthalpy

are:

dH,_, ~ al,- l,k

aHi-,

(12b)

(17b)

‘_’ ar,

ahi+, Ci.Ml.k= -L,+ 1------hi+,; ii.I+l.k c.I,MZ,Ml= For stages

(13b)

1 and N, the equations

The derivatives

UW

As already

are simply:

L?1.MZ.k - RN.,,., = 1. of y, are given as follows : yij = K;lij;

ah x = bj,. K$ - YsjKZ d + Yij

K; = Kij/x Ki,li, m

c:H.

+

aKij ahi ahi dK, shown, the derivatives~~~~~and

, ’ ,

are found by analytical

"_Vij

(14b)

differentiation

dH,

I

of the equation

of

state. The following equation relates derivatives with respect to molar flow rates to derivatives with respect to mole fraction :

Ldb4=ah4 c ,d~ axj i xs axi alj (15b)

I, ’ Ik

where M is In I#Qor AH.

(19b)