Ind. Eng. Chem. Fundam. 1983, 22, 472-485
472
Material Stability of Multicomponent Mixtures and the Multiplicity of Solutions to Phase-Equilibrium Equations. 1. Nonreacting Mixtures Davld B. Van Dongen5 and Mlchael F. Doherty' Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 0 1003
James
R. Halghtt
Department of Food Engineering, University of Massachusetts, Amherst, Msssachusetts 0 1003
The theory of material stability of multicomponent nonreacting mixtures is critically discussed and carefully presented using vector-matrix notation. Some singular transformations introduced in standard treatments are removed and several new forms of the stability and spinodal criteria are developed. The stability of two-phase mutticomponent mixtures is also discussed and it is found that there can be nonunlqueness of the equilibrium point. Some implications for the design of multicomponent azeotropic distillation systems are mentioned.
1. Introduction
The problem of determining whether or not a homogeneous multicomponent mixture will spontaneously and irreversibly split into two or more separate and distinct phases has been an important and active area of research since the problem was first addressed by Gibbs (1876), Duhem (1893), and Jouguet (1921). Even today, the literature contains several attempts to clarify the solution to the problem, such as those by Caram and Scriven (1976), Gray (1981), and Heidemann (1978). The criteria for material stability are well-known and have been welldocumented in many texts, e.g., Prigogine and Defay (1967) and Model1 and Reid (1974); yet there remain several areas where the theoretical development needs clarification. The theory as applied to binary mixtures is not in doubt and may even be considered trivial, but the extension of the theories to multicomponent mixtures is nontrivial and lacks careful development. In some texts, the criteria are modified, af'ter the fact, in order to match the known results. The need still exists for a mathematically rigorous approach to stability so that the theories can be applied to more complicated systems involving chemical reactions, interfaces, etc. The purpose of this paper is to carefully derive multicomponent stability criteria and to compare the results with those already in the literature. Because of the careful development, several new facts about material stability and the spinodal curve have been found. A stability analysis is then applied to a two-phase multicomponent mixture in order to establish the stability of solutions to the equilibrium equations. The results show that multiple solutions can be found and that the material stability criteria can help to determine which of these solutions may be stable. Preliminary Definitions. In the body of this paper, the discussion will center around an examination of the total Gibbs free energy, G, of a closed system relative to allowable changes that may take place within the system under the assumption of constant temperature and pressure. Classical analysis states that equilibrium will be Hunt-Wesson Foods, Inc., 1645 W. Valencia Drive, Fullerton, CA 92634. $Union Carbide Corp., P.O. Box 670, Bound Brook, NJ 08805. 0196-4313/83/1022-0472$01.50/0
attained when the total Gibbs free energy of the system attains a global minimum. If an equilibrium state is to be calculated, then some tests must be performed to ensure that the calculations do indeed predict a state that is at a global minimum of G. These tests can be quite complicated and require a tremendous effort in trial and error searches. For that reason, simpler stability criteria have been derived which guarantee that the calculations predict a state that is at least a local minimum of G and hence stable with respect to small perturbations. The total Gibbs free energy can be written in several ways depending on which description is more convenient for the topic in question. It is assumed to be an additive function of the total Gibbs free energies for all phases present in the system, and the total Gibbs free energy for each phase is assumed to be a first-order homogeneous function of the mole numbers of each species in that phase. The temperature, T , and pressure, P, are considered to be parameters in the problem and are held constant. The state vedor n p = (nlP,...,np) is introduced as a convenient shorthand notation for a description of the composition and total size of phase p . For a system of ?r phases, each containing all c components, the total Gibbs free energy is written as T
G = CGp(np,TQ) p=l
(1)
where each GP can be written as
C
GP = CnFpF(XP,T,P)
(3)
GP = [Cnlp]gP(XP,T,P) = nPgP(XP,T,P)
(4)
i=l
C
i=l
where n p is the total number of moles in phase p , p,P is the chemical potential of species i in phase p , and g p is the molar Gibbs free energy of phase p ; both are functions of the state vector of c - 1mole fractions XP = ( x l P , ..., x C J , T , and P. The first-order homogeneity of GP is important since it is the basis for the Gibbs-Duhem equations, which can be written as 0 1983 American Chemical Society
Ind. Eng. Chem. Fundam., Vol. 22,
No. 4, 1983 473
The total Gibbs free energy of the original mixture before the phase split will be GO = ng(x,TP) (12) and the total Gibbs free energy after the split will be
or
G = nlg(x + tl, TQ)+ n"g(x + &I, TQ) (13) It is assumed that the molar Gibbs free energy is an analytic function defined for all values of x near the composition of the original mixture. This means that g(x + c) can be expanded in a Taylor series up to second-order terms g ( x + c) =
or
and
where These equations can also be written in the compact matrix form where
&' is a c
X
&.np = 0 c square matrix
(9)
In ordef for eq 9 to have a nontrivial so_lutionfor np,the matrix GP must be singular, that is det (GP) = 0. This fact is almost always overlooked in classical stability analysis (see section 5). 2. Material Stability of a Homogeneous Mixture Some of the subject matter at the beginning of this section is well-known but is included for the sake of completeness of the presentation and, more importantly, because it introduces some of the procedures and notation which will be needed later. Suppose that a mixture of composition x exists as a single homogeneous phase a t a given pressure and at a temperature just below some critical temperature. Our first problem is to find the criteria which determine whether or not this mixture will spontaneously split into two or more phases. These criteria are called the criteria of material stability, sometimes called diffusional stability or intrinsic stability. If it is imagined that the homogeneous phase does split and in so doing decreases the total Gibbs free energy of the closed system, then the laws of thermodynamics state that the phase split may occur in reality. Whether or not the new heterogeneous system actually forms cannot be determined by thermodynamics alone. On the other hand, if the original phase is stable, then the formation of a heterogeneous system will increase the total Gibbs free energy, and thermodynamics guarantees that this process cannot occur under any circumstances. Since the mixture exists in a state just below the critical temperature, then two phases (denoted phase I and phase 11)may form out of the original mixture, each of which has a composition only slightly different from that of the original mixture. Let the new phases differ by small amounts c1 and cu where c1 is a state vector of c - 1deviations cf = xf - xi. Let there be n1and nn total moles in each of these are related by the overall new phases where the n's and ~i)s and species mass balances i = 1, ..., c - 1 (11) nxi = nl(xi E:) nn(xi + $1);
+
+
n = n1
+ n11
are evaluated at the original conditions of x, T,and P. Substituting these expansions into the expression for G and making use of the mass balances yields the change in the total Gibbs free energy during the process as a quadratic form 6G = G - Go = f/z[n*+ (15) The zero- and first-order terms have dropped out because of the mass balance relations. The change, 6G,can also be written in a compact matrix form as
6G = l/z[nl+ (n1)2/n11]t1T.gc1 (16) The matrix g is a (c - 1)X (c - 1)square symmetric matrix consisting of the second partial derivatives of g with resped to the mole fractions
If 6G is positive the phase split will not occur since it increases the total Gibbs free energy of the system. In consequence, the phase is materially stable if and only if 6G > 0. A necessary and sufficient condition for 6G to be positive for all small fluctuations is that the matrix g be positive definite. This, in turn, can be tested by Sylvester's theorem (see Mirskey, 19611,which states that a necessary and sufficient condition for a matrix to be positive definite is that its symmetric part have strictly positive principal minors. Some useful corollaries of this theorem are: (1) the matrix must have strictly positive elements on the main diagonal; (2) the determinant of the matrix must be strictly positive; (3) the eigenvalues of the matrix must be strictly positive. The limit of material stability of a homogeneous mixture is defined as the locus of points at which g just ceases to be positive definite. This locus of points is called the spinodal curve and, following the arguments of Mirskey (1961), this occurs where the determinant of g itself becomes zero. Thus, the spinodal curve is defined by det(g) = 0. For a binary mixture, see Figure 1, the criterion of material stability reduces to ( d 2 g / d x I 2 ) T g> 0. A plot of
474
Ind. Eng. Chem. Fundam., Vol. 22, No. 4,
1983 const P
I stable region
/crztical
I
point
06
.
5.4
t.
'3.2
.
x2
,"
I
-
C4
02
3
methyl ucetate
Figure 1. Relationship between the molar Gibbs free energy, g, composition, q,the binodal curve, and the spinodal curve for a binary mixture with a miscibility gap.
1
I
06
38
1C vethanol
XI
Figure 3. Relationship between composition, the binodal curve, and the spinodal curve for a mixture of methanol, hexane, and methyl acetate at 15 " C (288.16K).
Tie lines between the binodal curves are found by solving the set of nonlinear equations
piI(d,TQ)= p?'(d',T,P)(i = 1, 2, 3)
3
L
0
I
I
1
0.2
0.4
The relationship between the binodal curve and the spinodal curve is shown in Figure 3 for this ternary mixture. The nonideal liquid phase was modeled by the 2-suffix Margules activity-coefficient equation. All relevant thermodynamic properties can be found in the paper by Doherty and Perkins (1979). As expected, ideal solutions can be shown to be materially stable for all compositions using the above criteria. The molar Gibbs free energy of an ideal mixture is defined as
I
00
36
1.0 C
Tethonol
xi
methyl ocetote
gid
-----I,
+ RT In x i )
Cxi(p: i=l
Figure 2. Spinodal curves for a mixture of methanol, hexane, and methyl acetate at different temperatures.
[
("> (*,)
3x1' T , P , x , T,P
(9") axiax,
T,P
I
(5) T,P,x
=O
(18)
2 for a mixture of hexane, methanol, and methyl acetate at 1 atm pressure and at several different temperatures.
(20)
which produces the g matrix (eq 21), where zcis the de-
.. 1
+1
g vs. x1 will reveal whether or not the mixture exhibits a miscibility gap. If the curve is convex, the mixture is stable; if the curve is concave, the mixture is unstable and will never exist as a homogeneous phase. On the same plot, the binodal curve can be located by drawing a series of common tangents on the g - x1 curves. The spinodal curve is the locus of inflection points on the g - x1 curves and the region in between the two sets of points is called the metastable region. Homogeneous liquids can exist in this region but only under carefully controlled conditions. In ternary mixtures, the spinodal curve is found by solving eq 18. The solution of eq 18 is shown in Figure
det
(19)
gid
XC + 1 . x,
"
1
. . 1
XC xc- 1
I
t 1
pendent mole fraction. This matrix can be broken into separate parts (eq 22) where D is a diagonal matrix and
Ind. Eng. Chem. Fundam., Vol. 22,
Q is a matrix filled with 1’s. Because of the special form of this matrix, the determinant of gidis given by c-1 c-
c-1
det(gid) = (RT/Xc)c-l[~Xc/Xi+ r=l
c (dxc/xi)i+k]
(23)
&=I r = l
The determinants of the principal minors of gidare given by similar expressions with appropriately reduced indices. Since the mole fractions are always positive, the determinants of gidand its principal minora are all positive; hence g“‘ is a positive definite matrix and the ideal solution is materially stable for all values of x . An alternative but equivalent definition of the spinodal curve has been derived by Doherty and Kim (1981) which makes use of the expressions for the chemical potentials, pi, as functions of the mole fractions. The authors show that the spinodal curve can be found by solving the equation det(p) = 0 where p is a (c - 1) X (c - 1) square nonsymmetric matrix
No. 4, 1983 475
Therefore, along the spinodal curve, the determinants of g, p, and G all vanish and any one can be used as a definition of the curve. The relationship between matrices g and G can be strengthened. We will now prove that they are both positive definite for materially stable mixtures. This is a useful result which we will use repeatedly. Therefore, our objective is to prove that G is positive definite if and only if g is positive definite. To prove this, the molar Gibbs free energy is written explicitly in terms of the mole numbers, ni, and by using the notation C
n = Eni i=l
we can write G as G = ng(nl/n, $/n, ..., ne1/n). BY application of the chain rule for differentiation
and by further differentiation Still another equivalent definition of the spinodal curve can be found in terms of the total Gibbs free energy as a function of the mole numbers of all species in each phase. This approach is based on the definition of the chemical potential (understood here to be applied to each phase separately)
-L[( k=l
(5) dXjaXk
hence
TP,~!
This form may seem cumbersome but each term of G can be written in a more compact form involving the matrix g as a bilinear form
G, = l / n ( x - ei)T.g(x - ej)
A new matrix, G, can be defined as the submatrix of fi involving only the first c - 1 rows and columns and which is related to the p matrix by expression 27. This matrix - zXkPik
[‘;
P21 -
Pi2
- XxkPik
P22 - X X k P 2 k 2xkC(2k X X k P c - i , k Pc-1,2 - z X k P c - i , k
v P c - i ,~
Pi,c-1 - X X k P i k P2,c-i - zxkP2k Pc-i,c-1 - z X k P c - i , k
can be written in a more compact form as 1 1 G = ;(p - p-X) = -[p.(I - X)] n
3
(27)
L.c-l
: XC-1
. ...
I
..‘xc-lJ
The determinant of G is simply the product of the determinants of U and (I - X). According to Ferrers (1855), the determinant of (I - X) is simply x,; hence
(
det(G) = ;)‘-‘xc
detW
G = l / n ( X - I)T.g(X - I) To prove that G is positive definite, we need to show that the quadratic form cT.G.e is strictly positive for all c # 0. Thus 1 eT.G*e= -tT*(X - I)T*g(X- I)*t n 1 = -[(X - I).c]T.g[(X - I).e] n
(28)
where X is a matrix containing identical column vectors, each being the composition vector x
I:
where ei is a (c - 1)-dimensional vector with unity in the ith position and zero elsewhere. The matrix G can be written still more compactly as
where u = (X - 1)s~.Thus the quadratic form is strictly positive so long as g is positive definite and provided u # 0. If e # 0, then u can be zero only if (X - I) is singular, and this will occur if 1 is an eigenvalue of the matrix X. The matrix X has one nonzero eigenvalue equal to (1 - x,) and (c - 2) zero eigenvalues. Thus, the spectrum of X does not include unity. To show that the implication goes both ways, we merely invert the problem g = n[(X - I)-’lT.G-(X - I)-]
and so g is positive d e f i t e so long as G is positive definite. The matrix (X - I) is invertible since we have already shown that ita determinant is equal to x,, which is nonzero, and the proof is complete.
476
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983
3. A Sufficiency Condition for Stability Frequently, thermodynamic data for mixtures are represented by the excess molar quantities which can be readily used to estimate the stability of a multicomponent mixture. For example, the molar Gibbs free energy can be written as where gidis the molar Gibbs free energy of an ideal mixture and gEis the excess molar Gibbs free energy. Thus C
gid=
Cxilc~~(TP) + RT in xi]
i=l
(32)
6G =
n'g'(x + €1)
+ nIIgII(x + €11) + nIIIgIII(x + 211). - ng(x) (37)
and
+ ti') + n"(xi + E!)' + I~.'"(x;+ €2") = nxi
n'(xi
(38)
n = .I + + $1 As before, the values of g', glI, and g"' are replaced by a suitable Taylor series expansion which is simplified by making use of the mass-balance and equilibriwm constraints to give expression 39, where the vedors and matrix
i=l
where w/' is the chemical potential of species i in the pure state at temperature T and pressure P and yi is the activity coefficient of species i in the mixture. This definition can be inserted into the matrix g of eq 17 where
1 2
6G = -
(33)
and then eq 16 can be written as
have been conveniently partitioned into (c - 1)-dimensional units. The matrix g is the same as that defined earlier. The quantities N' and N" are merely used for shorthand notation
(34) As shown previously, gidis a positive definite matrix. Therefore it is sufficient to show that if the matrix gEis positive definite, then the mixture will be stable. The simplest form of nonideality in liquid mixtures is the two-suffix Margulea activity-coefficientmodel, which gives for ternary mixtures (Doherty and Perkins, 1978)
The above notation can be further simplified to read 1
6G = -rT.G1-r 2
(41)
We can gain more useful information by rearranging the quadratic form in order to isolate the matrix g. Let
Gz =
(42)
LI-GI
where The matrix gE will be positive definite if
(A13
+ A23 - 2 & G G )
A13
0
A23
0 for all possible variations and hence the two-phase azeotropic equilibrium will be stable. However, we also obtain the strange result that 6G can be positive even if one of the matrices G' or GI' fails to be positive definite, providing the composite matrix (G' + Gn) is positive definite. That is, the two-phase azeotropic equilibrium might be stable if one of the phases is materially unstable providing the other phase is so strongly stable that the composite matrix (GI + G") is positive definite. This odd result does not occur for nonazeotropic systems because of the form of eq 81. Further research in this area would be quite worthwhile, especially if accompanied with some testing of specific models for twophase azeotropic mixtures. 7. An Application of Stability Theory. Nonuniqueness of Solutions to Phase-Equilibrium Equations As we have shown in previous sections, the stability of an equilibrium two-phase binary system can be determined from the material stability of each individual phase. To begin, it will be useful to examine the connection between an equilibrium calculation, that is, the solution of the set of equations c~/[TP,n;/(n;
+ ni)] = ~iv[TQinlV/(nlV + n2')l /. + n; = Ni (i = 1, 2)
(88)
(that is, given T, P, N1,and N,, solve for the equilibrium values of nll, n i , nlv,and n2")and the corresponding geo-
480
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983 vapor
GL=O-
0
O K - - - -
\
-
a
x , Y,
f 1
Figure 5. Vaporjliquid equilibrium phase diagram for a binary mixture with no miscibility gap.
Figure 6. Cross-sectionalview of the G’ surface of a binary mixture along a line of constant composition.
metric representation of the equilibrium point on a total G surface. Recall that a stable equilibrium occurs at any point in the (n?, n i ) space of independent mole numbers which corresponds to a minimum in the total Gibbs free energy
G = G’(nll, ni) + Gv(nlv,nZv)
(89)
We have omitted T and P from the argument lists of G1 and Gv since they are to be understood as parameters in the problem. As will be seen, the application of the previously derived stability criteria allows us to conclude that multiple solutions to the equilibrium equations can be obtained under special circumstances. To begin the analysis, consider a simple binary liquidvapor system as shown by the T - x - y diagram in Figure 5. Let the system be closed so that the overall composition
-
0
x:
Figure 7. Cross-sectionalview of the G1surface of a binary mixture along a line of constant molar mass; special case where molar mass is unity.
is constant and also allow the system to have a constant pressure. Then, depending on the temperature of the system, we will find all vapor above temperature P, all liquid below temperature P , and liquid-vapor mixtures at temperatures in between. At fixed temperature there is a unique solution for the state of the system in terms of compositions of both phases and the relative proportions of liquid and vapor. Correspondingly, if we were to plot the total Gibbs free energy, G, as a function of n: and n i where
Since G’ is a fiist-order homogeneous function, then G‘ can also be written as G’ = (n; n2)g’(xl) = nTkl(z1) (93) A slice of the G’surface along a line of constant composition gives a straight line through the origin with slope g’(xl) as shown in Figure 6. If we were to slice the G’ surface along the line
G = G1(n,’, ni) + Gv(nlv,n2’)
nT1=I
nll
(94)
or
+ nlv = N ,
nll
+
nZ1 nZv= N z
(91)
then we should find only one minimum point on the G surface. In order to visualize the G surface, we first need to examine the G‘ and G’ surfaces individually and then observe how the two surfaces are added. In the case of interest, the liquid and vapor phases do not exhibit any miscibility gaps so a plot of g’ (or g’) vs. x (or y) for any given temperature will be everywhere convex. Let us examine in detail the shape of the G1 surface. Recall that G’ is given by GI = n:(plol + RT In x1 RT In yl/ + nJ(pzo’+ RT In xz + RT In y2) (92)
+
+
+ni = 1
then the result would be a graph of g’(xl) as shown in Figure 7. Together, these slices enable us to visualize the G1surface as shown in Figure 8. Since the vapor phase is assumed to be ideal, the Gv(nlv,nZv)surface should look much like that of Figure 8 and is calculated from the expression G’ = nlv(pIov+ RT In yl]+ nZV(p2”” + RT In y z ) where nlv, n2’,yl, and yzcan all be expressed in terms of the independent mole numbers nll and nJ using the material balances in eq 91. To add the two surfaces properly, we must set the origin of the G’ surfact at the point (Nl, N2)in the (n:, n2’>plane and have the axes travel parallel to those of the G’ surface but in opposite directions. This arrangement is required by the component mass balances. The result will be the
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983
t
T t
/J origin, G L = 0
GL
vopor
-
0
2
line of constant composition, x:
481
1
X I . Y1
Figure 10. Vapor/liquid equilibrium phase diagram for a binary mixture with a liquid-miscibility gap.
Figure 8. Schematic view of the G1surface showing the relationship between the cross-sectional views.
0
02
04
x,
06
08
hexane
10 methanol
Figure 11. Extension of the vapor/liquid equilibrium phase diagram into the liquid-miscibility gap.
Figure 9. Schematic construction of the G surface by adding the GI and Gv surfaces.
total Gibbs free energy as shown in Figure 9. At the absolute minimum of G, the values of n: and n: allow us to calculate the compositions of both phases and the relative proportions of the two phases, just as we could do from the T - x - y diagram. This minimum in G corresponds to a stable equilibrium point which implies that both phases are materially stable
Now let us examine a more complicated case where there exists a heterogeneous azeotrope. The T - x - y diagram for this system is shown schematically in Figure 10. By making use of this figure, we can completely determine the structure of the system knowing the temperature and overall composition. Above the azeotropic temperature, TA”, we should expect all liquid, all vapor, or liquid-vapor mixtures depending on the overall composition and tem-
perature. Below the azeotropic temperature, there will appear either homogeneous or heterogeneous liquid mixtures depending on the temperature and overall composition. However, below the azeotropic temperature, the nonlinearity of the equilibrium equations results in multiple equilibrium solutions as will be shown below. If we were to ignore the heterogeneity of the liquid phase and assume that the highly nonideal liquid is homogeneous for all compositions, then a careful and extensive set of solutions to the equilibrium equations yi =
PFt(T )
p x i r i
(i = 1, 2)
(96)
would lead to a T - x - y diagram as shown in Figure 11. If we extend the vapor-liquid envelopes on the T - x - y diagram down into the region of liquid heterogeneity, a very rich structure appears. In Figure 11, this is shown for the system methanol/hexane a t 1 a$m total pressure. The vapor phase is assumed to be ideal and the liquidphase activity coefficients were modeled using the twosuffix Margules equation. Also included in Figure 11 are the spinodal and binodal curves demarcating the boundaries of the unstable and metastable regions for the liquid. Before continuing, we must retreat slightly and examine the G1surface for such a nonideal liquid mixture. The G’
482
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983
t
liquid phase 1
atm
1
GL
GL
rnetas7able'
Figure 14. Schematic diagram of the G1surface of a partially miscible binary mixture. Figure 12. Detailed view of the GI surface for a mixture of methanol (1) and hexane (2) at 322.5 K.
for the mixture and has elements (dG/dn;])~p,~b (i = 1, 2, ..., c). It is well known that
This gradient system of differential equations (97) exhibits the following properties. (i) Singular points (i.e., steady-state points or equilibrium points) of eq 97 in the phase plane (n:, nJ) (or state space n for multicomponent mixtures) correspond exactly to stationary points in the G(n l) surface. (ii) The singular points of eq 97 are either stable nodes, unstable nodes, or saddles (i.e., the eigenvalues of the Jacobian matrix are always real). (iii) Minima in G(nI) give rise to stable nodes in the state space of eq 97; maxima give rise to unstable nodes; saddles give rise to saddle points. (iv) The function G(n') decreases along every trajectory, i.e.
'
/' "r'
'
Figure 13. Another view of the GI surface of Figure 12.
surface for the same mixture at 322.5 K is shown in Figures 12 and 13. The regions of instability and metastability are found by locating the inflection points of the gl(xl) plot and the points of common tangency as shown in Figure 1. These concepts are represented schematically in Figure 14. Addition of the two surfaces can lead to multiple local minima in the G surface but these are difficult to see graphically as the minima are relatively shallow compared to the surrounding terrain. It is easier to visualize the G surface and locate its minima by transforming the probem into a gradient system of ordinary differential equations as follows. Construct the set of nonlinear ordinary differential equations dn' _ - -VG(n')
(97) dr where VG(n l) is the gradient vector of the Gibbs surface
-
= VG(nl) dn'/dt = V G ( d ) (-VG(n')) = -(VG(n') VG(nl)) < 0
-
(99)
Therefore, trajectories of eq 97 approach either a minimum in G(n') or strike the boundary of the state space. In view of these properties, the state space of eq 97, which for our example is simply the phase plane (nll, nJ), gives a convenient representation of the G(n') surface. By property (ii) above, stable two-phase equilibrium points occur only at stable nodes in the phase plane. Before proceeding, we might mention that eq 97 have some physical significance, although this is incidental to their use in this work. Imagine a closed liquid-vapor system having arbitrary compositions in either phase and such that the mass balances in eq 91 are satisfied. If the system is not at a stable equilibrium, mass will transfer between the liquid and vapor phases in a direction of decreasing chemical potential so that the total Gibbs free energy decreases. The mass-exchange process continues until the total Gibbs free energy attains a minimum (at least locally). If both phases are completely mixed and the entire system is held at a constant temperature and
Ind. Eng. Chem. Fundam., Vol. 22, No. 4, 1983
483
Z 0
02
C4
n:
C5
08
10
Figure 15. Trajectories of eq 97 for a mixture of methanol (1)and hexane (2) at a temperature above the dew point.
pressure, then the number of moles in the liquid phase will change according to the equation dn)/dt = kA(p7 - ~ 2 )
(100)
where 12 is a mass-transfer coefficient, A is the instantaneous interfacial area between liquid and vapor phases, and the driving force for mass transfer is represented by the difference in chemical potentials for species i in the vapor and liquid phases. The chemical potentials at a fiied temperature and pressure are only functions of the composition of the respective phases and thus indirectly functions of n/. The quantity kA may not be constant but it will at least be always positive and so we may transform the system of eq 100 into a set of equations dni _ (i = 1, ...)c ) (101) d r p: - p; Since
-
0
x1. Y,
Figure l6. Locations of metastable vapor/liquid equilibria for a binary mixture with a miscibility gap; tie lines connect z and y, and r ’ and y’ at temperature T. solnodal
10
08
06 n:
04
-
02
(dG/dn,l)Tp,nL= p/ - Pi Then eq 101 can be written dnl/dr = -VG(n’) which is identical with eq 97. For our example, a mixture of hexane and methanol was chosen with an overall composition of 2 = 0.5 as shown in Figure 11. Above the temperature of about 324 K, the system should be all vapor and the minimum in the total Gibbs free energy should occur at the point nI1 = 0 and n i = 0. Solving eq 97 for a variety of initial conditions produces the phase-plane representation of the G(n l) surface shown in Figure 15. In this figure T > 324 K and all initial conditions give rise to trajectories which end at the origin. Now consider a temperature just slightly below the azeotropic temperature (see Figures ll and 16). We can use the extrapolated liquid-vapor envelopes to predict two solutions to the equilibrium equations. These correspond to the points x - y and x ’ - y’ in Figure 16. Indeed, if we examine the total G surface (see Figure 17), we find two local minima (stable nodes); one corresponding to each of the equilibrium compositions from Figure 16. Each of these points exists within the metastable region of the liquid composition simplex and as such the previously derived stability analysis shows that each of these equilibrium points is stable. If an experimental setup could be constructed so as to handle the liquid phase with extreme care in order to prevent heterogeneous nucleation, then it should be possible to experimentally verify the claim of multiple solutions to the equilibrium equations.
0 0
02
04
n:
06
38
10
Figure 17. Trajectories of eq 97 for the case shown in Figure 16 with overall composition 2 = 0.5: two local minima are found: (0) equilibrium point.
z
I
0
x1, Y1
10
Figure 18. Locations of metastable and unstable vapor/liquid equilibria for a binary mixture with a miscibility gap; tie lines connect r and y, x’and y’, #”and y”, and r’”and y”’at temperature T.
Because of the structure of the T - x - y diagram in Figure 11,we can find still more solutions to the equilibrium equations at lower temperatures. For example, if the temperature drops such that the overall system is within the swallowtail section of the T - x - y diagram, then two more solutions to the equilibrium equations appear. These are shown schematically in Figure 18 and the corre-
484
Ind. Eng. Chem. Fundam., Vol. 22,
No. 4, 1983 spinodal
0
32
C4
n)
06
08
tc
Figure 19. Trajectoriea of eq 97 for the case shown in Figure 18 with overall compoeition 2 = 0.5; two local minima, two saddle points,and equilibrium point. one constrained minimum are found: (0)
sponding trajectories to the gradient system are shown in Figure 19. The two new equilibrium points, x ” - y”and x ”’-y ”’,appear on the phase plane as saddle points. They are true equilibrium points insofar as they satisfy the equilibrium equations: p/ = p ; (i = 1 , 2 , ...,c). However, they are clearly unstable since they correspond to liquids which are materially unstable. Some of the trajectories in Figure 19 do not converge on any of the local stationary minima; rather, they converge on the point (Nl,N2). This means that one potential minimum point is a single homogeneous liquid phase. However, we should not be fooled into believing that this point is stable with respect to all possible fluctuations just because it is a local constrained minimum in G for the liquid-vapor mixture. If we apply the material stability test to this liquid phase we find that G’ (n: = Nl, n: = N2) is not positive definite, which conforms with the position of this liquid phase inside the spinodal region on Figure 18. If we replace the vapor phase with a second liquid phase and repeat our two-phase calculations with appropriate changes in the expression for the total Gibbs free energy of the system, we will now find that the (n?, n;) phase plane displays a unique stable node interior to the boundaries. The values of G at this node is lower than any of the G values for the nodes in Figure 19 and thus we conclude that at the given conditions of T , P, N,, and N 2 the system exhibits a stable two-phase liquid-liquid equilibrium point. The above analysis and discussion brings up a very important fact. The point (Nl,N 2 ) is a local constrained minimum in G for the liquid-vapor mixture and trajectories are locally attracted to this point. However, trajectories are repelled away from this same point if the two-phase system is changed to a liquid-liquid system. Therefore, the character of the phase plane and the interpretation of it depends critically on the types of phases present and great care should be taken to ensure that spurious phase planes are not accepted as true representations of the systems’s behavior. Finally, at sufficiently low temperatures (i.e., below the swallow tail region in Figure 18)the phase plane exhibits a minimum at the point ( N l ,N2)as shown in Figure 20. This point is also unstable with respect to liquid-liquid phase splitting and the above discussion about the same point in Figure 19 applies here also. It is worth making the final note that at fixed P, N,, and N2,the “bifurcation” temperature at which these liquid-
0
02
04
nk
06
08
10
Figure 20. Trajectories of eq 97 for the mixture shown in Figure 18 at a temperature below the minimum temperature on the bubble point curve; one constrained minimum is found.
vapor phase planes start to display a constrained minimum at (Nl,N2)is the temperature a t which the overall composition line intersects the bubble point curve in Figure 18. 8. Some Brief Remarks on Heterogeneous Azeotropic Distillation Recent studies have shown that multiple steady-state solutions can be found for the composition profile in heterogeneous azeotropic distillation columns (Magnussen et al., 1979; Prokopakis and Seider, 1980; Propopakis and Seider 1981). (Even though multiple column profiles have been found in simulations, their existence in real industrial azeotropic distillation processes has not been confirmed in the open literature.) At present this phenomenon is quite poorly understood, but it seems clear from the material presented in section 7 that one possible cause of this phenomenon is the presence of multiple metastable solutions to the phase equilibrium equations in the vicinity of heterogeneous azeotropes. We do not claim that multiple metastable equilibria are the only possible cause for such multiple column profiles, but we merely point out that the phase equilibrium calculations are likely to be a very sensitive component of the heterogeneous distillation problem. We are presently pursuing this in more detail and will report our results at a later date. Conclusions We have presented a careful and reasonably complete discussion on the stability of homogeneous and multiphase mixtures of many components. The results obtained are shown schematically in Figure 4. In section 7 of the paper, liquid-vapor phase-equilibrium calculations are formulated as a gradient system of ordinary differential equations. The trajectories in the phase plane are special projections of the total Gibbs free energy surface for the mixture with a one to one correspondence between extrema in the G surface and singular points in the phase plane. For nonideal mixtures, multiple equilibrium solutions (i.e., multiple steady states in the phase plane) have been found. The asymptotic stability of each steady state in the phase plane is directly related to the material stability of each phase in the mixture. Some possible implications of our findings for heterogeneous azeotropic distillations are briefly mentioned. Acknowledgment The authors wish to thank the donors of the Petroleum Research Fund administered by the American Chemical
Ind. Eng. Chem. Fundam., Vol. 22,
Society for providing financial support, and the University of Massachusetts Computing Center for providing computer time and services. We are also indebted to R. Willey of the University of Massachusetts for providing Figures 2 and 3 in the text. Nomenclature A = interfacial area Ai, = interaction coefficients in two-suffix Margules activity coefficient model c = number of components D = matrix filled with 1’s ei = vector with 1 in ith position, 0 elsewhere G = total Gibbs free energy Gif = (a2GP/dnipdnf)Tg,,t GP = (c - 1) X (c - 1) matrix of elements Gi,P G j = linear transform or exploded form of g GP = c x c matrix of elements Gif g = molar Gibbs free energy gij = (a2g/axiaxj)Tg,.f g = (c - 1) X (c - 1)matrix of elements gi, H, Hi = linear transform or exploded form of GP I = identity matrix i , j , k, 1 = component indices k = mass-transfer coefficient Li = linear operators Ni = total moles of i in the system Np = ratios of mole numbers in different phases NP = vector of elements nip/np n = total moles in the system nip = total moles of i in phase p np = vector of elements nip P = pressure Pimt = saturation vapor pressure of pure i Q = diagonal matrix of elements l / x i R = gas constant T = temperature t = time X = matrix with x as column vectors x = mole fraction in liquid phase x = vector of elements x i y = mole fraction in vapor phase Z = overall composition for the system Greek Letters a = vector of deviations
activity coefficient 6 = difference operator aij = Kronecker delta function tip = deviations xip - x i y =
c
q
No. 4, 1983 485
= vector of deviations
= vector of deviations
B = vector of deviations p = chemical potential 7r
= total number of phases
u = vector of deviations 7 T
p
= time = vector of deviations
= matrix of partial derivatives, (dpCL,/dxj)T,p,x,
Superscripts I, 11, p , 1, v = designations of a particular phase id = refers to ideal solutions E = refers to excess molar quantities -1 = inverse operator T = transpose operator ’ = part of a larger vector 0 = reference state Subscripts i , j , k, 1 = component indices 0 = initial conditions
Literature Cited Caram. H. S.; Scrlven. L. E. Chem. Eng. Sci. 1976, 3 1 , 163-168. Doherty, M. F.; Kim, Y. S. Chem. f n g . Sci. 1981, 36. 1260-1263. Doherty, M. F.; Perklns, J. D. Chem. Eng. Sci. 1978, 33. 589-578. Doherty, M. F.; Perkins, J. D. Inst. Chem. Eng. Symp. Ser. 1979. No. 56, 4,2121-4.2148. Duhem, P. “Travaux et MBmoires des FacultBs de Lilie”, 1893; Vol. I11 B. Ferrers, M. Quart. J . Meih. 1855. 1 . 364. Glbbs, J. W. Trans. Connect/cut Acad. 1876, 3 , 108-248; 1878, 3 , 343-524. Gray, P. Chem. Phys. Len. 1981, 81. 457-480. Heldemann, R. A. Chem. Eng. Sci. 1878, 33. 1517-1528. Jouguet. E. J . .€cole Po/ytech. (Park),2nd Ser. 1821, 21, 81. Magnussen, T.; Mlchelson. M. L.; Fredenslund, A. Inst. Chem. Eng. Symp. Ser. 1979, No. 56, 4.2f1-4.2119. Mlrskey, L. “An Introduction to Linear Algebra”; Oxford University Press, London, 1961; pp 400-421. Modell, M.; Reid, R. C. “Thermodynamicsand its Applications”;PrenticaHall, Englewood Cliffs, NJ, 1974; pp 188-230. Prigoine, I.; Defay. R. “Chemlcal Thermodynamics”; Longmans, London, -1967; pp 205-228. Prokopakls, G. J.; Seider, W. D. Presented at the 73rd Annual Meeting of AIChE, Chlcago, 1980. Prokopakis, G. J.; SeMer, W. D. Presented at the 74th Annual Meeting of AIChE, New Orleans, 1981.
Received for review April 2, 1982 Revised manuscript received February 2, 1983 Accepted July 11, 1983
Presented at the 1982 Annual Meeting of AIChE, Los Angeles, CA, Nov 17, 1982.
