Fluid Phase Equilibria 158–160 Ž1999.

617–626

State function based flash specifications

)
Michael L. Michelsen
Institut for Kemiteknik, Bygning 229, DTU, DK 2800 Lyngby, Denmark
Received 6 April 1998; accepted 30 September 1998

Abstract

A variety of flash specifications of practical importance can be formulated as minimization of a thermody-

namic state function. It is well known that the solution for the PT-flash yields the global minimum of the
mixture Gibbs energy, but in addition specifications of PH, PS, TV, UV or SV all permit selections of
thermodynamic state functions for which a global minimum must be located. Two important advantages are
obtained with the minimization based approach. First, the desired solution is known to be unique, and second,
stability analysis can be used to verify its correctness and to determine the number of equilibrium phases. For
the PT-flash the solution can be determined by unconstrained minimization, whereas the remaining specifica-
tions are accompanied by one or two nonlinear constraints and thus less straightforward to attack. We present
here two approaches for dealing with such specifications. The first is a nested optimization approach where T
andror P are the dependent variables in the outer loop, and where a PT-flash is solved in the inner loop. The
essential advantage of this formulation is that it is very easy to implement but the drawback is the additional
cost of the nested loops. The second approach is based on a modified objective function in which all constraints
are removed, but where a saddle point rather than a minimum must be located. The resulting equations are
solved by a global Newton’s method, and it is shown that a common Jacobian matrix can be used with all the
specifications given above. q 1999 Elsevier Science B.V. All rights reserved.

Keywords: Method of calculation; State function; Equation of state; Gibbs energy

1. Introduction

A number of phase equilibrium calculations share the important characteristic that the solution
corresponds to the state that yields the global minimum of a thermodynamic state function. As a
consequence, except for a few degenerate cases like, e.g., azeotropes where the phase amount can be
indeterminate, these problems have a unique correct solution. The classical, and most frequent

)
Tel.: q45-42-88-3288; fax: q45-42-88-2258

0378-3812r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved.
PII: S 0 3 7 8 - 3 8 1 2 Ž 9 9 . 0 0 0 9 2 - 8
618 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626

example is the isothermal flash Ž P and T specified. where Gibbs energy minimization coupled with
tangent plane stability analysis has been widely used to determine the number and composition of the
equilibrium phases. Other traditional examples are the isenthalpic Ž PH . and the isentropic Ž PS . flash,
and in addition specification of volume and temperature Žstorage vessels and pipeline shutdown. and
specification of volume and internal energy Ž unsteady state operations. are becoming of increasing
importance.
In the present work we investigate whether the principles used for the isothermal flash can be
utilized for the remaining state function based specifications.

2. State functions

In Table 1 below are listed the state functions to be minimized with the corresponding flash
specifications.
We consider in all cases a feed of 1 mol of composition z. For the two-phase PT-flash the desired
molar flows z and l are found as the solution to:

min G Ž T , P ,z,l . Ž1.

subject to the constraints

T s T spec , P s P spec , z q l s z Ž2.

With T and P given, elimination of, e.g., the liquid phase flows reduces Eq. Ž 1. to the unconstrained
minimization problem:

min G Ž z, z y z . Ž3.
in which we only require that all compositions are nonnegative. Furthermore, it is only necessary to
determine a local minimum of G. Afterwards, stability analysis by means of the tangent plane
criterion enables us to verify whether the global minimum has been located. If not, new phases have
to be introduced, and the local unconstrained minimization is repeated.
The determination of a local, unconstrained minimum can be performed efficiently and routinely
w1x with a very high reliability by requiring that each iterative step must reduce the objective function.
The more difficult step is that of performing the stability analysis w2,3x, which requires a global
minimization. For two-phase calculations this is unproblematic w3,4x and for multi-phase calculations
heuristic methods have been applied with success, and more rigorous approaches are being investi-
gated w5,6x.

Table 1
State function to be minimized for a given specification
Specification P,T P, H P,S T,V U,V S,V
State function G yS H A yS U
 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626 619

The situation is however not as straightforward for the remaining specifications. The minimization
formulation for, e.g., the isentropic flash is:
min H Ž T , P ,z,l . Ž4.
subject to the constraints:
P s P spec , l q z s z , S Ž T , P ,z,l . y S spec s 0 Ž5.
Elimination of the linear material balance constraint as for the PT-flash results in:
min H Ž T , P spec ,z, z y z . Ž6.
subject to:
S Ž T , P spec ,z, z y z . y S spec s 0 Ž7.
The constraint of specified entropy is however nonlinear in the independent variables, and the
constraint cannot be eliminated explicitly. In principle, one could use the same procedure as for the
PT-flash, but the complicating factor in the constrained minimization is that the individual steps can
no longer be tested for a reduction in the objective function. Therefore, we do not have the same
guarantee of convergence as for the PT-flash.

3. Modified objective functions

All the above specifications can be transformed to a form that formally eliminates the constraints
and enables us to use the ‘natural’ variables T, P and the phase flows as independent variables.
Continuing with the isentropic flash as our example we introduce the function Q defined by:
Q Ž T , P spec ,z, z y z . s G q TS spec Ž8.
where the independent variables are T and z. The gradient vector is:
EQ EQ
s yS q S spec , s mzi y m il , i s 1,2, . . . ,C Ž9.
ET Ezi
and we observe that at the desired solution the gradient of Q with respect to all the independent
variables equals zero. The solution is thus a stationary point of Q, but unfortunately not a minimum
since Ž E 2 Q . rŽ ET 2 . s yŽ Cp .rŽ T . - 0, and the Hessian matrix is therefore not positive definite, as
required for a minimum. Similar Q-functions, all with the Gibbs energy as the ‘core function’, can be
formulated for the other specifications, as summarized below in Table 2.

Table 2
Q-functions for state function based specifications
Specification Q-function
P, H Ž Gy H spec .r T
P,S GqTS spec
T,V Gy PV spec
U,V Ž GyU spec y PV spec .r T
S,V GqTS spec y PV spec
620 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626

It is worthwhile mentioning that it may be possible to define modified Q-functions for which the
desired solution is the unconstrained minimum. One example is the isentropic flash with:
2
Q mod s G q TS spec q a Ž S y S spec . Ž 10.
where a is a constant, chosen sufficiently large. Procedures based on this choice, together with a
similar formulation for the isenthalpic flash, have been investigated by Michelsen w7x but are not
entirely unproblematic.

4. Maximizing Q

The desired stationary point of the above Q-functions are all saddlepoints of the Q-surface, having
positive curvature in the composition directions and negative curvature in the T- andror P-direction.
This opens up the possibility to nest a PT-flash calculation Ži.e., an unconstrained optimization in the
composition variables. with a maximization with respect to the remaining variables. For the isentropic
flash we get:
max Q s Ž Gmin q TS spec . Ž 11.
with respect to T, where Gmin is the solution to min GŽ T, P,z, z y z . at the current temperature. In a
similar manner all the other specifications can be solved by maximizing Q, combined with an inner
loop minimization of the Gibbs energy at the current temperature and pressure, i.e., an isothermal
flash.
If only a single phase is known to form at the specified condition the inner loop is unnecessary, and
the outer loop calculation is therefore essentially equivalent to a single phase equilibrium calculation.
The flash performed in the inner loop is readily differentiated with respect to temperature and
pressure to provide the gradient and the Hessian required for the determination of the maximum. For
the isentropic flash we obtain:
EQ EGmin
s q S spec s ySmin q S spec Ž 12.
ET ET
where Smin is available from the current phase distribution. The second derivative is found from:

E2 Q ESmin Cp,min
2
sy sy Ž 13.
ET ET T
but its evaluation requires the derivatives of the flash solution since:
EH EHl EHz C EHl El i EHz Ezi
Cp,min s ž / ž / ž /
ET P,z
s
ET P ,l
q
ET P ,z
q Ý
is1
ž El i ET
q
Ezi ET /
C EHz EHl Ezi
s Cp,l q Cp ,z q Ý
is1
ž Ezi
y
El i / ET
Ž 14.
 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626 621

where Cp,l are Cp,z the total heat capacities for the equilibrium phases. It is the last contribution that
requires information about the change in the equilibrium composition with temperature.
The main advantage of the nested loop approach is that it enables us to handle a variety of
specifications in a fairly simple manner, provided an efficient and reliable PT-flash is available. The
outer loop iteration involves only a single or two independent variables, and reliability is assured by
the use of a maximization approach. For the isentropic and the isenthalpic flash, a simpler and
probably widely used procedure is to utilize the PT-flash in an inner loop and to treat the temperature
determination in the outer loop as a root-finding, rather than as an optimization, problem. The
additional cost-free information from the underlying optimization formulation, however, provides for
faster convergence as well as increased robustness. Finally, the extension to multiphase calculations,
possibly also involving multiple solid phases, only requires the corresponding multiphase abilities of
the PT-flash.

5. A Newton approach

An obvious drawback of the nested loop approach is its lack of efficiency as compared with a
simultaneous convergence of all independent variables. When good initial estimates are available a
Newton-based approach may therefore be preferable. For the two-phase equilibrium calculation we
can derive a general formulation that is capable of handling all the above specifications with a
common Jacobian matrix. As independent variables are used the vapour phase flows, z, lnT and ln P
Žprovided one or both are not specified.. The Newton iteration uses the general formulation:

M gT gP g
Dz

 g TT
g PT
ET T
ET P
ET P
EPP 0 0  0
DlnT q r T s 0
Dln P rP
Ž 15.

where only the elements r T and rP depend on the specification. The deviation vector is given by:

g i s ln yi q ln w i Ž y,T , P . y ln x i y ln w i Ž x ,T , P . , i s 1,2, . . . ,C Ž 16.

with the two variable elements being:

Specification rT rP
Ž T, P . – –
Ž P, H . Ž H spec y H .rŽ RT . –
Ž P,S . Ž S spec y S .rR –
Ž T,V . – P Ž V y V spec .rŽ RT .
Ž V,U . Ž U spec q PV spec y H .rŽ RT . P Ž V y V spec .rŽ RT .
Ž V,S . Ž S spec y S .rR P Ž V y V spec .rŽ RT .
622 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626

and the elements of the symmetric Jacobian matrix are:

E gi
Mi j s , i s 1,2, . . . ,C, j s 1,2, . . . ,C,
EÕj

Eln w i Ž y,T , P . Eln w i Ž x ,T , P .

g T ,i s T ž ET
y
ET / , i s 1,2, . . . ,C

Eln w i Ž y,T , P . Eln w i Ž x ,T , P .

g P ,i s P ž EP
y
EP / , i s 1,2, . . . ,C, Ž 17.

Cp P EV P 2 EV
ET T s y , ET P s , EPP s
R R ET RT EP
where Cp and V are the heat capacity and the volume of the combined phases.

6. Helmholtz function based expressions

The Ž T,V . -flash can also be solved as an unconstrained minimization of the Helmholtz energy if
we abandon the traditional calculation of thermodynamic properties at specified temperature, pressure
and composition in favour of a calculation at specified temperature, Õolume and composition Žwhich
is in fact the ‘natural’ choice of independent variables for a pressure-explicit EOS. . This formulation
is advantageous in itself for critical point calculations and in connection with chemical models or
models with density-dependent mixing rules. It is interesting to note that in addition to the VT-flash
also the PT-flash can be carried out using unconstrained minimization with ‘ volume-based’ thermo-
dynamics. The objective function

Q s A q VP spec Ž 18.
where z, Vz and Vl are chosen as the independent variables, satisfies at its stationary point the
equilibrium conditions
EQ EQ EQ
s mzi y m il , s yPz q P spec , s yPl q P spec Ž 19.
Ezi EVz EVl

and in addition, the second derivatives have the correct sign,

E2 Q EPz E2 Q EPl
2
sy ) 0, 2
sy )0 Ž 20 .
EVz EVz EVl EVl

The ‘ volume-based’ PT-flash appears first to have been proposed and investigated by Nagarajan et al.
w8x.
 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626 623

Table 3
Q-functions using temperature and volume as the independent variables
Specification Q-function
P,T AqVP spec
P, H 1r T Ž AqVP spec y H spec .
P,S AqTS spec qVP spec
T,V A
U,V 1r T Ž AyU spec .
S,V AqTS spec

Q-functions for the remaining specifications with the Helmholtz energy as the ‘core function’, i.e.,
based on temperature and volume as the independent variables, are listed in Table 3. Only for the
specifications Ž T,V . and Ž T, P . does the solution correspond to a minimum of Q.

7. Solution strategy

When no initial estimates are available the Wilson K-factor expressions can be used to generate
approximate values of phase properties. Assuming the vapour phase to be ideal, we obtain
Pc i Tc i
ln w iz s 0, ln w il s ln K iWilson s ln ž / q 5.373 Ž 1 q v i . 1 y ž / Ž 21.
P T
Combined with an expression for the ideal gas heat capacity as a function of temperature this enables
us to calculate all relevant thermodynamic properties. Obviously, the inner loop calculation of phase
compositions reduces to solving the Rachford–Rice equation for the vapour fraction.
The approximations obtained by means of the Wilson approximation is normally adequate when
the vapour fraction at the solution is not very small. The most obvious defect is the severe
misrepresentation of liquid volumes, which are evaluated to identically zero. Its main advantage is
that it enables a simple and safe determination of an initial approximation to the solution.
Our recommended approach for solving the flash equations can be summarized as follows.
Ž1. Converge the equations, using the Wilson K-factor based relations. The resulting estimate of
temperature, pressure phase fraction and phase composition cannot be expected to be very accurate
but will often prove ‘reasonable’.
Ž2. Continue with a few Ž5–10. steps of a partial Newton’s method based on the equations derived
above, but neglecting the composition derivatives of the fugacity coefficients. The partial Newton’s
method is likely to underestimate the magnitude of the corrections, and this may be advantageous at
an early stage, in particular when the initial estimates are inaccurate. When composition derivatives
are neglected, the structure of M enables us to reduce the effective size of the set of equations to 3 or
less. For moderately nonideal mixtures this successive substitution like approach may well converge
the set of equations in a modest number of iterations.
Ž3. Attempt to converge the set of equations using the full Newton’s method Ž including explicitly
the composition dependence of the fugacity coefficients in the Jacobian. .
Ž4. In case of a failure in step 2 or step 3, e.g., lack of convergence, excessive corrections in
temperature or pressure or removal of a phase, switch to the Ž safe. nested loop approach.
624 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626

Table 4
Iteration history for four test cases
Solution T r P T and P from T and P after five steps Full Newton
ŽKrMPa. Wilson approximation with partial Newton iterations
240.00r4.052 220.6r4.01 240.00r4.052 None
210.00r4.052 203.1r4.832 210.01r4.052 2
220.00r5.065 211.08r6.047 220.01r5.065 3
200.00r5.065 212.9r7.412 200.57r5.080 4

Ž5. Check the final result for stability, and introduce new phases if required.
Step 4 is rarely necessary, but the availability of a slower but reliable procedure enables us to take
full advantage of the efficient but more risky Newton-based approach.

8. Examples

The approach outlined above is tested with the energy–volume specification on a seven-component
natural gas mixture containing 94.3% methane, 2.7% ethane, 0.74% propane, 0.49% n-butane, 0.27%

Fig. 1. Q-function for UV-flash at equilibrium point Ž240 K, 4.05 MPa..

 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626 625

Fig. 2. Q-function for UV-flash at equilibrium point Ž200 K, 5.07 MPa..

n-pentane, 0.1% n-hexane and 1.4% nitrogen, using the SRK-equation. The mixture critical point is
203.1 K, 5.89 MPa. Test phases were generated by carrying out PT-flash calculations at a number of
points and subsequently using the values of U and V from the converged solution as specifications for
the UV-flash. The results are shown in Table 4, where T and P from the Wilson approximation and
after five steps with the partial Newton method are given, in addition to the subsequent number of full
Newton steps required. None of these examples required the use of step 4..
For two of the test cases, contour plots of the Q-function, together with the phase envelope for the
seven-component mixture are shown in Figs. 1 and 2. It is worthwhile to notice that the contours of
constant Q intersect the phase boundary without a change in slope. The contours in the Ž 240 K, 4.052
MPa. case of Fig. 1, which is characterized by a high vapour fraction, are very regular, indicating that
the maximum is easily located, whereas the contours of Fig. 2, close to the bubble line, change more
abruptly. For this case, the vapour fraction changes rapidly with the temperature.

9. Conclusion

A formal framework has been set up that enables us to formulate a variety of phase equilibrium
problems as unconstrained maximization problems by means of a PT-flash calculation in an inner
loop. As an efficient alternative, a general form of Newton’s method can be used, with the
maximization as a backup for difficult cases.
626 M.L. Michelsenr Fluid Phase Equilibria 158–160 (1999) 617–626

An interesting variant, not yet investigated, is a formulation in terms of volume and temperature
and with the Helmholtz energy as the core function.

10. Nomenclature

A Helmholtz energy
ET T , ET P , EPP Jacobian matrix elements
G Gibbs energy
gT , gP Jacobian matrix elements
H Enthalpy
Ki Equilibrium factor, component i
l Liquid moles vector
M Matrix of composition derivatives
P Pressure
Q Objective function
R The gas constant
r T , rP Deviation vector elements
S Entropy
T Temperature
U Internal energy
V Volume
z Vapour moles vector
x, y, z Mole fraction vectors
Greek letters
mi Chemical potential, component i
wi Fugacity coefficient of component i

References

w1x R. Fletcher, Practical Methods of Optimization: Vol. 1. Unconstrained Optimization, Wiley, New York, 1980.
w2x L.E. Baker, A. Pierce, K.D. Luks, Soc. Petrol. Eng. J. 22 Ž1982. 731–742.
w3x M.L. Michelsen, Fluid Phase Equilibria 9 Ž1982. 1–19.
w4x M.L. Michelsen, Fluid Phase Equilibria 9 Ž1982. 20–39.
w5x C.M. McDonald, C.A. Floudas, AIChE J. 41 Ž1995. 1798–1814.
w6x J.Z. Hua, J.F. Brennecke, M.A. Stadtherr, Fluid Phase Equilibria 116 Ž1996. 52–59.
w7x M.L. Michelsen, Fluid Phase Equilibria 33 Ž1987. 13–27.
w8x N.R. Nagarajan, A.S. Cullick, A. Griewank, Fluid Phase Equilibria 62 Ž1991. 191–210.