THREE PHASE FLASH CALCULATIONS FREE ENERGY MINIMISATION
USING
M. E. SOARES* and A. G. MEDINA Centro de Engenharia Quimica, Faculdade de Engenharia, Porto, Portugal
C. MCDERMOTT and N. ASHTON Department of Chemical Engineering, University of Birmingham,Birmingham, England (Receiued 25 May 1981: accepted 21 August 1981)) Abslnct-The distillation of systems exhibiting partial miscibility over part of the composition range is growing in interest and the modeling of three phase distillation columns is highly desirable. This paper presents work which has been produced as part of a project directed towards such objective. An algorithm for flash calculations in which partially miscible systems are involved is described. The algorithm is based on the principle of minimisation of the system Gibbs free energy which is becoming increasingly popular allowing as it does for the simultaneous satisfaction of the constraints of equilibria and mass balance requirements. The minimisation procedure is based on a Newton-Baphson algorithm and so in order to avoid trivial solutions it is necessary to have an adequate choice of initial conditions. In the present paper this is taken into due account and a policy for selecting initial compositions and flowrates is presented. Numerical examoles are aresented for the ternary systems ethanol/ethylacetate/water and methanollethylacetate/water. . INTRODUCTION
The importance of liquid-liquid and liquid-liquid-vapour flash calculations in process design has lead to the appearance of several papers during recent years. The method of minimisation of the system Gibbs free energy has proved to be increasingly popular allowing as it does for the simultaneous satisfaction of the constraints of chemical and physical equilibria and mass balance requirements. Among the early suggestions for direct application of free energy minimisation was the description by Dluzniewski et al. (1973) of a comprehensive programme which included the possibility of multiple liquid phases, chemical reaction equilibria and the use of Lagrange multipliers to handle the constraints. A similar procedure qescribed by George et al. (1976) removes some of the difficulties caused by restricted ranges of values for some of the independent variables. The method involves a simple change of variables using an allocation function so that the constrained problem is reexpressed as an equivalent unconstrained problem. Essentially the same procedure is proposed by Heyen et ol. (1979); in this work an initial application of the simplex algorithm is included to establish ranges of values for the independent variables in order to preclude negativity problems in the presence of constraints arising from reaction equilibria. A common problem with the above methods is the necessity to specify in advance the total number of phases to be expected in the solution. This information may not be available in the case of partially miscible *Author to whom correspondence should be addressed.
systems where the feed composition may lie in regions of either one, two or even three liquid phases. At the specified conditions a vapour phase may or may not be present. This problem is discussed by Heidemann (1974) as applied to the minimisation of free energy of liquidliquid systems in the absence of chemical reaction. Here the free energy minimum is located by the integration of a set of differential equations whose right hand sides are the derivatives of free energy with respect to phase component flowrates. In this form the possibility of negative phase flowrates is allowed and these are interpreted as absence of the component from the phase or absence of the phase itself. Heidemann also adopts a change of variable to work with the logarithm of the flowrate when it is known that the flowrate can not become negative. This is also coupled with a choice of independent variable to ensure that the phase flowrate decreases as the calculation progresses. By this means mole fractions greater than unity can be avoided. A further possibility which arises with partially miscible systems is that a feed composition may lie in the metastable region as discussed by Prausnitz and Maurer (1979) and by McDermott er aI. (1980) among others. A feed with composition in this region may exist as a single phase or if subjected to sufficiently large disturbances will split to form two or more liquid phases. Although the possibility of phase splitting under these circumstances may be recognised from a plot of free energy for binary systems as described by Heidemann (1974) no explicit account is taken of this possibility in the minimisation procedure. Unfortunately no local stability test, capable of detecting metastability as opposed to total stability is currently 521
M. E.
522
SOARES
available. However McDermott et al. (1980) have shown that by suitable choice of initial guessed compositions, the trivial solution is more easily avoided. It must be recognised that, in the case of a metastable feed, the trivial solution is physically realizable and so forms a valid solution for the minimisation procedure. However in engineering applications it is likely to be the solution with liquid phase splitting which is of interest. The calculation method proposed in the present paper for three phase flash is an extension of the method of McDermott et al. (1980) to include the possibility of a vapour phase. CALCULATION
The method presented flash calculations is based tion of Gibbs free energy For a system with m phases, the free energy is
where ni = total molar flowrate of component j; P,’ = vapour pressure of component j; P ; total pressure;yii = activity coefficient of component j in phase i; xii= mole fraction of component j in phase i; pF=chemical potential of component j (standard state). Equation (2) can be differentiated with respect to the molar flowrate of each component in each phase taking into account that the total number of moles of each component is specified. Considering the vapour phase molar flowrates as the dependent variables we have aGIRT ani,
i=, ,z = & - ga
(3)
where g,, = In P;Y+~~; g3i = In Px,.
(4)
The minimisation procedure is based on a NewtonRaphson algorithm in which the Marquardt method (Marquardt, 1963) is used to control the correction vector (HfLI)an=-g where
an, correction
vector;
matrix; I, identity matrix; g, first derivatives vector; A, variable that controls the correction vector. The use of the above method requires the evaluation of second derivatives. These can be written as:
where
METHOD
in this paper for three phase on the principle of minimisa(McDermott et nl., 1980). components and y different given by
where G = Gibbs free energy of the system; ng = molar flowrate of component j in phase i: pij =chemical potential of component j in phase i. For a system with two liquid phases and a vapour phase, considered as ideal, the above expression reduces to:
w
et al.
(5) H,
second
derivatives
and
H3j,x, = L (&Z R,t x3i
l);&
= l(i = j); S,, = 0
(i# j).
t-1
The use of the Marquardt method to control the magnitude of the correction vector calculated from eqn (5) also has the effect of removing problems caused by any ill-conditioning of the matrix, H, of second derivatives from eqn (6). Such singularities can be expected when two liquid phases tend towards the same composition leading to a linear dependence of the matrix H. It has proved useful to make also use of the gradient method to avoid too large composition changes above specified limits, which could lead to trivial solutions. The use of eqn (5) for isothermal flash calculations requires an adequate choice of initial values of the independent variables. In previous work involving liquid-liquid equilibria calculations it was found that, in order to avoid trivial solutions, the initial guess should overestimate the separation to he obtained. This could he achieved considering each phase i as being almost pure component i with small specified amounts of the other components. However if a vapour phase is present a different approach can he adopted and this is exemplified in Fig. 1 which refers to a type I system. The choice of initial conditions starts with an estimation of the composition and flowrate of the vapour phase (point V’) hased on calculated values of boiling and dew point temperatures of the feed at the specified pressure (TB and TD). The calculation of boiling point temperatures for mixtures either in the two phase or in the single phase regions requires the assumption of two phases which overestimate the separation to be obtained and the choice of an initial guess for temperature. Phase compositions are then calculated using the liquid-liquid equilibrium calculation method described by McDermott et nl. (1980). Such method assumes a number of liquid phases equal to the number of components. This number is eventually reduced as phases with equal compositions or zero Rowrates will appear in the final solution. For any of the liquid phases obtained mole fractions in the vapour phase can he evaluated and if they do not add up to unity the Newton method is used to generate a new value of temperature. Convergence is obtained in a small number of iterations.
Three phase flash calculations using free energy minimisation NUMERKAL
R
L /”
C
ltESUL.TS
An algorithm based on the method described in the
“‘,
LA
523
-‘+‘-
7
_ -
-
- 4 0
Fii. 1. Definition of initial composition and flowrates: F, feed composition: L, global liquid composition; L,,L*, estimation of liquid compositions: V. vapour composition at the boiling point temperature of the feed; V’. estimation of vapour composition.
A linear interpolation with respect to temperature is carried out knowing that for T = Te the vapour flowrate is equal to zero and that for T = To the vapout flowrate is equal to the feed flowrate. Also for T = TB the vapour composition is the one obtained for the bubble point calculation (point V) and for T= TD the vapour composition is the feed composition (point F). This provides the first estimate of the composition and flowrate of the vapour phase at the specified system temperature T. A simple mass balance then leads to the corresponding overall liquid phase composition and flowrate (point L). This may then be split into two liquid phases (points L, and L) which overestimate the degree of separation to be expected. The above procedure allows the determination of an initial estimate of the number of moles of the different components in the different phases which is close to the solution. Thus a small number of iterations is normally required. The method described does not require the total number of liquid phases to be specified. If the liquid obtained as a result is in the single phase region the two liquid phases will have the same composition; alternatively one of the phases will be associated with a zero flowrate.
Table 2. Anloine constants (log P’ =
previous paragraph was developed; numerical results based on the use of the NRTL equations (Renon and Prausnitz, 1968) were obtained for the ternary systems ethanol/ethylacetate/water and methanollethylacetate/water studied by Zandijcke et al. (1974). For such systems Zandijcke et al. (1974) propose the set of parameters presented in Table 1. Vaponr pressures for the pure components can be evaluated in terms of the Antoine equation (Antoine constants are presented in Table 2). Calculated values of equilibrium compositions at boiling point temperatures are presented in Fig. 2 for the system ethanollethylacetate/water. Initial compositions in a two phase region and in a single phase region are
considered. In Table 3 calculated results are compared with experimental determinations of liquid-liquid-vapour equilibrium compositions as published by Zandijcke et al. (1974). Although the differences obtained are small (Table 4) there are difficulties involved in the determination of suitable sets of parameters for the correlation of liquid-liquid-vapour equilibrium and in fact for the simultaneous correlation of liquid-vaponr and liquid-liquid data. Flash calculations are exemplified in Fig. 3 which refers to a ternary mixture with 10.7% ethanol, 30.1% ethylacetate and 59.2% water (molar percentages) flashing at 71.3o”C and one atmosphere. For that mixture calculated values of bubble and dew point temperatures are 71.23”C and 8666°C respectively. Two different initial estimate and calculated values of phase composition and flowrates are presented in Table 5. It is apparent that the initial guesses based on the interpolation method previously described are close enough to the solution for convergence to occur in a small number of iterations. For the same mixture tlash calculations were performed at different temperatures. Results are presented in Fig. 4 where the total number of moles in the vapour
A - B/( C+
T); P*-mm
Hg; T-Y)
524
M. E. SOAIWSet al. ETHANOL
.1
.5
.6
.?
.9
ER
ID E~~ACETATE
Fig. 2. Phase compositions at boiling point temperatures. (Ternary system ethanoliethylacetatelwater; 0, liquid; A, vapour.
P = 1 atm).
ETHANOL
IMUiO n
_ 0
.i
.b
.I
.c
bATW
Fig.
3.
.s
.6
\
.1
.(I
.s
ID
ETWLACETATE
Phase compusitiuns fur different iterations. (Feed composition: ethanol-10.7%, ethylacetate-30.1%6, water-59.2%; T = 71.3VC;P = 1atm).
phase is shown as a function of temperature. This figure clearly indicates that three phases only exist in a narrow temperature range close to the boiling point. Flash calculations can be used to generate equilibrium diagrams for partially miscible systems. One such diagram (system methanol/ethylacetate/watcr: P=
1 atm: T = 70.6O”C) is presented in Fig. 5. Single phase, two phase and three phase regions can be identified. The shape of the ternary diagram depends strongly on the behavior of the three component binary systems. For binary systems it is possible to represent graphically (Fig. 6-8) the variation of Gibbs free energy of mixing
Three phase flash calculations using free energy minimisation
Fig. 4. The influence of flash temperature on the amount of the vapour phase. (Feed composition: ethanol-l&7%, ethylacetate30.1%. water-59.2%; P = 1atm).
.2
.‘
WTER
ETHYIACETATE
Fig. 5. Phase diagram for the ternary system methanoliethylacetate/water (T = 70.6CtY; P = 1 atm).
with composition considering the system vapour. For a binary system it can be written 1!974) that: AGM ~=“‘W’+(l-x,)ln~~
as liquid or (Heidemann,
(7)
where AG”, Gibbs free energy of mixing; fi, fugacity of component i. Equation (7) implies that a standard state with pure components at unity fugacity at the system temperature is assumed and that fugacities for the liquid and vapour
M. E. SONUS et
526
al.
Table 4. Root mean square deviation of calculated and experimental values for the ternary system ethanol (I)/ethylacetate
(2)lwater
0)
Table 5. Initial and final values for three phase flash (feed composition: ethanol 10.7%; ethylacetate 30.1%; water 59.2%; temperature = 71.3o”C)
.z
.A
mdc
.6
fmctim
.I
of LmybDaate
Fig. 6. Gibbs free energy vs composition for the binary system ethylacetatelwater (T = 70.6OC; P = I atm).
0
.I
.G
I?&
‘&mm eLlanu”
Fig. 7. Gibbs free energy YS composition for the binary system methanollethylacetate (T = 70.6O”C;P = I atm).
Three phase flash calculations using free energy minimisation IIETHANOC
WATER
ETHYLACETATE
Fig. 9. Phase diagram for the ternary system methanol/ethylacetate/water (T = 73.o”C: P = I atm).
mde Fractionof
m&hand
Fig. 8. Gibbs free energy vs composition for the binary system methanol/water (T = 70.6O”C:P = 1 atm).
phases
are calculated
as:
For each binary system the relative position of liquid and vapour curves clearly indicates whether one, two or three phases will exist. For ethylacetatelwater (Fig. 6) it is clear that at 70.6O”C and 1 atmosphere no vapour phase exists and that for ethylacetate mole fractions between 0.05 and 0.75 phase splitting occurs. For the binary methanol/ethaylacetate at the same temperature and pressure (Fig. 7) it is possible to draw a common tangent line to the liquid and vapour curves. This indicates that there is a one-liquid phase region (methanol mole fractions smaller than O.lO), a vapour phase region (methanol mole fractions greater than 0.26) and a two phase region. A similar reasoning can be applied to the binary methanol/water (Fig. 8), which presents a two phase region (liquid and vapour) for methanol mole fractions between 0.61 and 0.85. A ternary diagram for the same system at a different temperature (T = 73.0-C) is presented in Fig. 9. CONCLUSIONS
In the present work a new method for three phase isothermal flash calculations based on the principle of minimisation of the system Gibbs free energy was The minimisation procedure is based on a presented. Newton-Raphson algorithm in which the Marquardt method is used to control the correction vector. A method for selecting initial values of phase compositions and phase flowrates was proposed. This involves an
estimation of bubble and dew point temperatures of the feed and linear interpolations with respect to temperature and flowrates. In order to avoid trivial solutions the degree of separation to he expected is overestimated. Numerical examples obtained for the ternary systems ethanol/ethylacetate/water and methanol/ethylacetate/ water were presented. Convergence in a small number of iterations was obtained for all the situations that were considered. Acknowledgement-The authors want to express their gratitude to the Scientific Affairs Division of Nato for the award of Nato Research Grant No. 1715. The support of FundaCBo Calouste Gulbenkian-Lisboa-Porb& is also gratefully acknowledged.
NOTATION
4B,C Antoine constants NRTL oarameters
“/ fuga& derivative of free energy i first vector of first derivatives of free energy
I3 Gibbs free energy AGM H I m n n
Gibbs free energy of mixing square matrix of second derivatives of free energy identity matrix number of components number of moles vector of independent number of moles total pressure p9 pure component vapour pressure number of phases ; gas constant T absolute temperature x mole fraction
Greek Symbols a NRTL parameter S Diiac function y activity coefficient A constant used in Marquardtmethod p chemical potential Subscripts i phase i j component j
M. E. SOARESel al.
528 s r
phase s phase t
Superscripts 0 standard state L liquid phase V vapour phase REFERENCES Dluznicwski J., Adler S. Chem. Engng Prog. 1973 69 (11) 79. George B., Brown L., Farmer C., Buthod P. and Manning F. Znd. Engng Chem. Proc. Des. Ihev. 1976 G(3) 372.
Heidemann R.,, A.1.C.h.E. I. 1974 200) 847. Heyen G., Kahtventzeff B., Germain A. and Lecorsais R. Proc., Computer Applications in Chemical Engineering, Montreux 1979. Marquardt D., J. Sot. Znd. Appl. Marhs 1963 11431. McDermott C.. Ashton N., Soares M. E. and Medina A. G., Proc. Kemtek ‘S-5th Int. Gong. in Scandinavia on Chemical Engineering, Copenhagen, 1980. Prausnitz J. and Maurcr G., Znt. Symp. Distillation. Institution of Chemical Engineers, London 1979. Zandijcke F. and Verhoeye L., J. Appl. Chem. Eiotechnol. 1974 24 709.