SciELO - Scientific Electronic Library Online

vol.16 issue3Inulinase from Kluyveromyces marxianus: culture medium composition and enzyme extractionUsing hybrid neural models to describe supercritical fluid extraction processes author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Brazilian Journal of Chemical Engineering

Print version ISSN 0104-6632On-line version ISSN 1678-4383

Braz. J. Chem. Eng. vol.16 n.3 São Paulo Sept. 1999 

Phase rule calculations and the thermodynamics of reactive systems under chemical equilibrium


Curso de Pós-Graduação em Tecnologia de Processos Químicos e Bioquímicos (TPQB), Departamento de Engenharia Química, Escola de Química, Bloco E, Centro de Tecnologia, Universidade Federal do Rio de Janeiro, CEP 21949-900, Phone: (021) 590-3192, Fax: (021) 590-4991, Rio de Janeiro - RJ, Brazil, E-mail:, or


(Received: November 12, 1998; Accepted: May 12, 1999)



Abstract – In this paper, we examine the resolution of some phase rule problems within the context of multiple chemical equilibrium reactions, using cubic equations of state and an activity coefficient model. Bubble and dew reactive surfaces, reactive azeotropic loci and reactive critical loci are generated and presented in graphical form. Also isobaric bubble and dew reactive enthalpy loci, which may be useful in the modeling of reactive distillation operations, are depicted. All the formalism here employed is developed within the coordinate transformation of Ung and Doherty, which is appropriate for equilibrium reactive or multireactive systems. The major contribution of this work is the determination of critical loci for reactive or multireactive equilibrium systems. Since it is known that for some class of chemical reactions the kinetics and product distribution exhibit high sensitivity to pressure near criticality, the present study may be useful as a predicting tool in these cases if the chemical equilibrium condition is not too far from the real phenomenon.
Keywords: Vapor-liquid equilibrium, critical point calculations, transformed coordinates, reactive azeotropes.




The problem of vapor-liquid equilibrium with chemical reactions using activity coefficient models has frequently been addressed in the chemical engineering literature in recent years (Barbosa and Doherty, 1988a, 1988b, 1988c; Ung and Doherty, 1995a, 1995b; Thiel et al., 1997). This problem can be resolved by the specification (determined by the phase rule) of transformed coordinates (instead of molar fractions), proposed by Barbosa and Doherty (1988b) and Ung and Doherty (1995a). However, in this transformed context, approaches using cubic equations of state (CEOS) have not yet been utilized. Also the examination of reactive or multireactive equilibrium loci of first, second or higher order is a subject not commonly found in these previous works.

In this work, we present a multi-purpose phase rule specification formalism, appropriate for generation of equilibrium multireactive loci of first order (i.e., which are dependent on first derivatives of free energy models). This formalism is used with CEOS models (i.e., Helmholtz free energy models) and is adequate for the treatment of multireactive cases. The formalism is based on the Ung-Doherty coordinate transformation and is used in the generation of bubble, dew and azeotropic reactive loci for two systems of importance in the chemical engineering literature.

Critical point calculations received considerable attention in the chemical engineering literature in the 80’s (de Medeiros, 1990; Heidemann and Khalil, 1980), but the influence of equilibrium chemical reactions in the critical phenomenon was not considered. Recently, some experimental works on thermal decomposition of hydrocarbons (Yu and Eser, 1995, 1997a, 1997b) have pointed out that the kinetics and product distribution are very sensitive to pressure if the process occurs in the vicinity of a critical state of the mixture. If, additionally, one may assume that the chemical equilibrium condition is not an unrealistic description of the phenomenon, product distribution and reaction yields may be predicted with a suitable algorithm for equilibrium calculations of reactive criticality. In the present work, the problem of calculation of the critical point in reactive systems is solved within the Ung-Doherty context of compositional coordinates. An algorithm suitable to multireactive criticality, predicted with CEOS models, is presented, and typical critical reactive loci are generated.

CEOS applications developed here concentrate on two CEOS models: the van der Waals (vdW) (Walas, 1985) and the Soave-Redlich-Kwong (SRK) (Soave, 1972) equations. Although systems studied here have a well-known non-ideal liquid phase, vdW and SRK models were used with binary parameters set to zero because (i) as far as we can see, there is no published matter concerning the CEOS binary parameters needed here; (ii) statistical estimation of these parameters, though straightforward with vapor-liquid equilibrium data, is beyond the scope of this paper; and (iii) the paper focuses on formalisms and algorithms for exploration of multireactive equilibrium loci with CEOS models and is not particularly intended to modify, enhance or tune such models aiming at accurate prediction of these loci in comparison with available experimental information. Concerning our reasons for using vdW and SRK models, it may be said that (i) the SRK equation is very popular in chemical engineering applications; (ii) in spite of not having found any relevant usage in engineering, the vdW equation is the simplest CEOS, and frequent pedagogic applications were found; and (iii) the free-energy context employed [AV(T, r)] is not common in the chemical engineering literature (in fact, we suspect it has not yet been used), so we think it may be interesting for someone else to deal with the AV(T, r) formalism constructed for two CEOS: the simplest and another with widespread usage.

Finally, we would like to emphasize that all formalisms and algorithms presented in this work are completely general in order to be used with multireactive equilibrium problems. Unfortunately, due to a lack of reliable and industrially relevant published data on multireactive distillation or multireactive vapor-liquid equilibrium, all systems illustrated here are monoreactive but originated in two industrial processes: the MTBE system and the ethyl acetate system. It seems pointless to work with completely artificial multireactive systems, made with anonymous species, in order to merely demonstrate a multireactive application. Thus, in the entire text, the word "reactive" has a multireactive connotation when applied to algorithms, equations, formalisms and thermodynamic states.

Throughout the paper some abbreviations will appear in order to try to enhance the conciseness of the text. So, reactive vapor-liquid equilibrium, reactive azeotrope and reactive critical point will be referred to as RVLE, RAz and RCP, respectively.



Reactive Bubble and Dew Loci Using Activity Coefficient Models

In order to introduce the basic Ung-Doherty formalism and the systems which will be used in the following sections, standard vapor-liquid calculations under chemical equilibrium are presented. Liquid phase is assumed to be non-ideal and is described by the Wilson model (Smith et. al, 1996):

; (1)

The vapor phase is assumed to be an ideal gas. Use of activity coefficient models for vapor-liquid equilibrium problems with chemical reactions was previously addressed by Barbosa and Doherty (1988a, 1988b, 1988c) and Ung and Doherty (1995a, 1995b).

In the classic form, the equations that describe the chemical and phase equilibrium of a vapor-liquid system with nc components and nr chemical reactions are

(nc equations) (2)

; (3)

(nr equations) (4)

where operation "ln" produces an application of the natural logarithm to all components of a vector, thereby generating a vector of identical size; is the vector of vapor-liquid equilibrium factors [k = k(x, T, P)]; K is the vector of chemical equilibrium constants [K = K(T)] and represents the matrix of stoichiometric coefficients (the size of vector 1 is nc x 1).

This problem has nc + nr + 2 equations (equations 2-4) and 2nc + 2 variables (x, y, T and P). The phase rule shows that the number of degrees of freedom is nc – nr. Considering the system pressure previously set, the result is nc – nr – 1 degrees of freedom. Thus, one can specify nc – nr 1 components of the composition vector of the liquid phase or of the vapor phase in order to explore reactive bubble and dew loci, respectively.

Ung and Doherty (1995a) proposed a set of transformed coordinates (X, Y) – henceforth referred as UD coordinates –, based on the vectors of mole fractions, as follows: considering a system with nc components and nr equilibrium chemical reactions, nr reference components, whose subset (or subvector) is denoted by the designator "k" are chosen; thus, the UD coordinates for liquid and vapor phases are given, respectively, by





size(X) = size(Y) = (nc-nr x 1);

size(Xk) = size(Yk) = size(xk) = size(yk) = (nr x 1)

Accordingly, vectors of mole fractions are partitioned as

; (7)

With equations 5 and 6, it is clear that


It is valid to state that the choice of reference components must be made in order to guarantee the existence of matrix . Thus, inert components can not be chosen as reference components because they will decrease the rank of matrix .

It can be shown that, as in the case of mole fractions, the sum of UD coordinates is also equal to one, or

and (9)

These relationships can be verified directly from equations 5 by premultiplication by 1t. Thus, the reactive transformation of Ung-Doherty produces nc – nr – 1 independent compositional coordinates. Therefore, the reactive vapor-liquid equilibrium can be solved after specification of the vector of the UD independent coordinates. The UD coordinates can also be expressed in terms of molar densities of liquid and vapor phases, as explained later in this work.

Phase rule considerations and general implications concerning UD coordinates are shown in Tables 1, 2 and 3 for the reactive systems studied in this work. In these tables, the constraints defining the UD compositional domain, and thus its shape, are also presented. The phase rule balance of degrees of freedom (DF) is determined from the classic form, DF = nc + 2 - np - nr -nac, where nc, np, nr and nac stand for the numbers of components, phases, independent chemical reactions and independent additional constraints, respectively (for RVLE np = 2, nac = 0; for RAz np = 2, nac = nac - nr - 1; and for RCP np = 1, nac = 2).

Data necessary for resolution of RVLE for MTBE - I system are presented in Tables 4 and 5. Vapor pressures are calculated by the Wagner equation (Smith et al., 1996). The chemical equilibrium constant is given below (Barbosa and Doherty, 1988a):


where E1 = -4254.05, E2 = 10.0982 and E3 = 0.2667. The standard state is the (l) state, i.e., pure liquid at 1 bar and temperature T. The Newton-Raphson method (Ortega and Rheinboldt, 1970) is used to solve the set of non-linear equations (equations 2-4) after inclusion of (nc - nr - 1) equations 5 for specification of the independent UD coordinates. With pressure fixed at 1 bar, the bubble RVLE locus and the dew RVLE locus are surfaces plotted against UD coordinates in Figure 1.


Table 1: Metyl-tert-butyl ether I (MTBE I) system

Components: isobutene (1), methanol (2), MTBE (3), n-butane (4)
Reaction: (1) + (2) Û (3)
Reference component: MTBE (3)
Liquid phase UD coordinates:
; ; ;
UD compositional domain: 0 £ X1, X2, X4 £ 1; X1X2X4 = 1
Domain shape: triangle
Liquid phase independent UD coordinates: X1, X2
Domain vertices:

Isobutene (x1 = 1)Þ X1 = 1;  X2 = 0

Methanol (x2 = 1)Þ X1 = 0;  X2 = 1

n-butane (x4 = 1)Þ X1 = 0;  X2 = 0

Degrees of freedom and state specification:

RVLE: DF = 3, (P, X1, X2) or (P, Y1, Y2)

RAz: DF = 1, (P)

RCP: DF = 2, (X1, X2)



Table 2: Metyl-tert-butyl ether II (MTBE II) system

Components: isobutene (1), methanol (2), MTBE (3)
Reaction: (1) + (2) Û (3)
Reference component: MTBE (3)
Liquid phase UD coordinates:
; ;
UD compositional domain: 0 £ X1, X2 £ 1; X1 + X2 = 1
Domain shape: line segment [0,1]
Liquid phase independent UD coordinates: X1
Domain vertices:

Isobutene (x1=1)Þ X1=1

Methanol (x2=1)Þ X1=0

Degrees of freedom and state specification:

RVLE: DF = 2, (P, X1) or (P, Y1)

Raz: DF = 1, (P)

RCP: DF = 1, (X1)



Table 3: Ethyl acetate (EA) system

Components: acetic acid (1), ethanol (2), ethyl acetate (3), water (4)
Reaction: (1) + (2) Û (3) + (4)
Reference component: ethyl-acetate (3)
Liquid phase UD coordinates:
X1 = (x
1 + x3); X2 = (x2 + x3); X3 = 0; X4 = x4 - x3
UD compositional domain: 0 £ X1, X2 £ 1; -1 £ X4 £ 1; X1 + X2 + X4 = 1
Domain shape: square
Liquid phase independent UD coordinates: X1, X2
Domain vertices:

Acetic acid (x1 = 1)Þ X1 = 1; X2 = 0

Ethanol (x2 = 1)Þ X1 = 0; X2 = 1

Ethyl acetate (x3 = 1)Þ X1 = 1; X2 = 1

Water (x4 = 1)Þ X1 = 0; X2 = 0

Degrees of freedom and state specification:

RVLE: , or

RAz: ,




Table 4: Data for the MTBE - I and MTBE - II systems


Isobutene (1)

Methanol (2)

MTBE (3)

n-butane (4)


























pvB (K)(b)





pvC (K)(b)






93.33 x 10-6

44.44 x 10-6

118.8 x 10-6

100.39x 10-6

A (J/gmol.K)(d)





B (J/gmol.K2)(d)

2.804 x 10-1

7.092 x 10-2

5.136 x 10-1


C (J/gmol.K3) (d)

-1.091x 10-4

2.587 x 10-5

-2.596 x 10-4


D (J/gmol.K4) (d)

9.098 x 10-9

-2.852x 10-8

4.303 x 10-8



-1.691 x 104

-2.013 x 105

-2.931x 105


Tc (K)(d)





Pc (bar)(d)





w (d)





(a) Aznar and Telles (1993); (b) Barbosa and Doherty (1988a);
(c) Ung and Doherty (1995a); (d) Reid et al. (1987)



Table 5: Binary parameters of the Wilson model for the MTBE - I and MTBE - II systems (Ung and Doherty, 1995a) (K)




























Figure 1: Equilibrium diagram for the MTBE - I system at 1 bar with the Wilson model.


Generation of Transformed Enthalpy Surfaces in Multireactive Systems

The UD coordinates can be used to obtain the reactive extension of the McCabe-Thiele approximation, as presented by Platt (1997). In reactive cases, the McCabe-Thiele equations for enthalpies of saturated phases are

; (11)

where X and Y correspond to the UD coordinates, defined by equations 5. Transformed molar enthalpies (in the Ung-Doherty context) are defined as follows (Platt, 1997):

; (12)

Thus, equations 11 result in:

; (13)

We call equations 13 the reactive extension of the McCabe-Thiele approximation. These linear and parallel forms are useful in reactive distillation columns design because they promote the decoupling of material and energy balances, allowing the determination of compositional profiles with only reactive mass balance equations around sections of the column (Platt, 1997). Thus, it is worth checking the suitability of these equations for some reactive systems.

At low pressures, and without dimerization reactions, the residual terms of vapor enthalpies can be removed, and the following expressions are obtained for vapor and liquid molar enthalpies at reactive saturation:



where Tdew and Tbubble are reactive dew and bubble saturation temperatures, obtained at the (P,Y) and (P,X) coordinates, respectively. Dividing equations 14 and 15 by the denominators of equations 12, transformed enthalpies can be rigorously calculated and compared with equations 13.

Algorithm 1 is used to generate the loci of transformed molar enthalpies of saturated phases for the MTBE - II and EA systems. Liquid phase is modeled with the Wilson equation and vapor phase is assumed to be an ideal gas. In equations 14 and 15 latent heats for pure components are estimated with the Pitzer method (Reid et al., 1987) (J/gmol), as follows:




Algorithm 1: Generation of dew and bubble transformed enthalpies for reactive systems

1. System is chosen at specified pressure;
2. X is specified (nc - nr - 1 specifications);
3. Problem of reactive bubble point is solved (equations 2-4), and x, y and T are obtained;
4. Bubble enthalpy is obtained with equation 15, and transformed bubble enthalpy is calculated;
5. Steps 2 - 4 are iterated for entire transformed compositional space;

6. Diagram X versus transformed bubble enthalpy is generated;
7. Steps 2 - 6 are repeated for the problem of dew reactive point calculation.


Specific heats at constant pressure for liquids are calculated by (Reid et al., 1987)


and ideal gas heat capacities at constant pressure are obtained by (Reid et al., 1987)

image191.gif (1034 bytes)(18)

Vapor pressures are predicted by the Antoine equation (Barbosa and Doherty, 1988a) for both systems. Coefficients for the Antoine equation, molar volumes of saturated liquid (at 298.15 K), coefficients for specific heat at constant pressure for ideal gas, enthalpy of formation (ideal gas at 298.15 K) and critical properties for the EA system are presented in Table 6. Binary parameters for the Wilson model are shown in Table 7. Data for chemical equilibrium constant (equations 10) are E1 = -450.382, E2 = -0.8135 and E3 = 0 (Barbosa and Doherty, 1988c). Data for the MTBE - II system was presented in the previous section (Tables 4 and 5).


Table 6: Data for the EA system


Acetic Acid (1)

Ethanol (2)

Ethyl Acetate (3)

Water (4)

pvA (K)(a)





pvB (K)(a)





pvC (K)(a)






57.54x 10-6

58.69x 10-6

98.49x 10-6

18.07x 10-6

A (J/gmol.K)(b)





B (J/gmol.K2)(b)

2.549x 10-1

2.141x 10-1

4.072x 10-1

1.924x 10-3

C (J/gmol.K3)(b)

-1.753x 10-4

-8.390x 10-5

-2.092x 10-4

1.055x 10-5

D (J/gmol.K4)(b)

4.949x 10-8

1.373x 10-9

2.855x 10-8

-3.596x 10-9

(J/gmol) (b)

-4.351x 105

-2.350x 105

-4.432x 105

-2.420x 105

Tc (K) (b)





Pc (bar) (b)





w (b)





(a) Barbosa and Doherty (1988a); (b) Reid et al. (1987)



Table 7: Binary parameters of the Wilson model for the EA system (Barbosa e Doherty, 1988a) (K)





























Figure 2: Bubble and dew transformed molar enthalpies (J/gmol) for the MTBE - II system at 1 bar.



Figure 3: Bubble and dew transformed molar enthalpies (J/gmol) for the EA system at 1 bar.


In Figure 2, it is clear that the curves of bubble and dew transformed enthalpies for the MTBE - II system are approximately linear and parallel, validating the reactive McCabe-Thiele approximations for this system. On the other hand, as can be seen in Figure 3, the reactive McCabe-Thiele approximation does not correctly represent the surfaces of transformed enthalpies for the EA system.



A Specification Formalism for the Determination of Multireactive Loci of First Order

RVLE and RAz calculations can be performed using conventional CEOS, with a phase rule formalism for the specification of functions of the UD coordinates. Results are illustrated for the vdW equation and the SRK equation.

The specification formalism is based on equations describing the problems of RVLE and RAz, in the following equivalent format:

(nc equations) (19)

(nr equations) (20)



; (23)



where X and Y now stand for the independent portion of the vector of the UD coordinates [size (nc – nr – 1 x 1)], Av refers to the Helmholtz free energy per unit of volume and the chemical potential vector is defined by


where the derivatives are assumed to be at a constant temperature. Equations 21 and 22 establish the Ung-Doherty compositional transformation from molar densities vectors () for liquid and vapor reactive phases (matrices and are defined in the following text; see also equations 30). In these equations, is the molar density of component "i", K is the vector of chemical equilibrium constants and CEOS represents the cubic equation of state applied in the context. Equations 24 and 25 represent the necessary specification equations (as dictated by the phase rule) for several problems, as shown in Table 8.

The UD coordinates can also be defined from molar density vectors. In order to present this new form, we define the following selection matrices:

; (27)

where denotes a matrix where all elements are zeros and refers to an identity matrix. The sizes of these submatrices are (nc – nr – 1 x nc – nr 1), (nr x nr), (nc – nr – 1 x nr), (nc – nr – 1 x 1) and (nr x 1). Thus, the sizes of matrices and are (nc – nr – 1 x nc) and (nr x nc), respectively. With equations 27, the molar densities vector and the stoichiometric matrix can be written in the following form:

; ; (28)

where the following partitions are defined:

; (29)

and the sizes of entities , , are (nc - nr - 1 x 1), (nr x 1), (nc - nr - 1 x nr) and (nr x nr), respectively.

Thus, with partitions 29 and equations 27 and 28, the independent UD coordinates (for the liquid and vapor phases) may be described by

; (30)

Table 8 shows the degrees of freedom analysis and the necessary specifications for the following phase rule problems: (i) bubble RVLE at fixed UD composition; (ii) bubble RVLE at fixed pressure; (iii) dew RVLE at fixed UD composition; (iv) dew RVLE at fixed pressure; and (v) RAz. In Table 8, the superscript "mov" refers to a "moving" specification and "fix" refers to a fixed specification.

Generation of RVLE and RAz Loci Using Cubic Equations of State

Table 9 presents formulae necessary to establish the system of equations 19 - 25 for determination of molar density vectors of vapor and liquid phases, transformed compositional vectors X and Y, and temperature. These equations are solved by the Newton-Raphson method. For the vdW equation, we used the Lorentz-Berthelot mixing rules (Walas, 1985), and for the SRK equation we used the same mixing rules with the Graboski and Daubert coefficients (1978). Other definitions follow.

Molar densities:

image232.gif (951 bytes)  (31)



where is the chemical potential of pure component "i" as an ideal gas at P = 1 bar and temperature T.

In CEOS context and according to equation 32, the reference (or standard) state for the formation of chemical species is the (g) standard state, i.e., pure ideal gas at 1 bar and temperature T. Since the expressions for presented by Barbosa and Doherty (1988a) use (l) standard state (pure liquid at 1 bar and temperature T), the following transformation is useful:

; (33)


where "0" refers to pure ideal gas at P = 1 bar and temperature T and "*" corresponds to pure liquid at P = 1 bar and temperature T. In the calculations performed in this paper, it is assumed that







where j is the reaction index.

As an illustration, results are obtained for the MTBE - I system. Figure 4 shows the equilibrium diagram at 1 bar (bubble and dew RVLE surfaces) for this system, with the SRK CEOS, qualitatively similar to that obtained by the Wilson model (Figure 1), as expected. Figures 8 and 9 show bubble RVLE surfaces at 1 bar and 10 bar for this system, using the vdW and SRK CEOS. In this case, Psat (T) values are calculated using the Wagner equation (Table 4).

Determination of RAz locus is presented for the EA system, using the SRK equation. The analysis of degrees of freedom for this problem is presented in Table 8. Necessary data are shown in Table 6.

In order to illustrate the determination of the reactive azeotrope path with pressure, the concept of reactive residue curve map (RRCM), explained by Barbosa and Doherty (1998b), is invoked. The equations that describe the dynamics of simple distillation (also called "open evaporation") of reactive mixtures under specified pressure are


where and refer to the independent UD coordinates for the liquid and vapor phases, respectively, and q is a dimensionless time. On the RRCM, a reactive residue curve (RRC) can be obtained integrating 40 from a given initial condition. At each point on an RCC, a bubble RVLE problem must be solved for (T, Y).

A reactive azeotrope is a singular point of equation 40, that is, an RVLE condition such that


The dynamic stability character of such a point can be inferred from the RRCM by the behavior of neighbour RRCs. With 5, it is easy to see that non-reactive azeotropes and pure component states also satisfy 41 and are singular points on the RRCM, as well. On the RRCM, an increase in temperature on an RCC is usually indicated with arrows.

Figures 5 and 6 present the RRCM for the EA system at 1 bar and 13 bar, with the SRK equation. In these figures one can see that a quaternary reactive azeotrope (a global minimum temperature reactive azeotrope, i.e., an unstable node), which appears at 1 bar, moves to the ethanol / water edge at 13 bar. Figure 7 presents the RCP surface for this system and the path of this quaternary RAz. All calculations are performed with the SRK equation. The generation of the RCP surface, shown in Figure 7, will be explained in the next section.


Table 8: Degrees of freedom and specifications for some problems on reactive vapor-liquid equilibrium with CEOS

Phase-rule problem

Bubble RVLE

Bubble RVLE




Degrees of freedom

nc nr

nc – nr

nc – nr

nc – nr



X (fixed), P (moving)

X (moving), P (fixed)

Y (fixed), P (moving)

Y (moving), P (fixed)

P (moving)



nc – nr – 1


nc – nr – 1




Table 9: Formulae for the vdW and SRK equations of state


vdW equation

SRK equation

classical form



mixing rules

; ;






Table 9 (cont.)



Figure 4: Equilibrium diagram for the MTBE - I system at 1 bar with the SRK equation.



Figure 5: RRCM for the EA system at 1 bar with the SRK equation.



Figure 6: RRCM for the EA system at 13 bar with the SRK equation.



Figure 7: Path of quaternary reactive azeotrope and RCP surface for the EA system with the SRK equation.



Determination of critical coordinates in reactive systems is approached within the Ung-Doherty context, and programs were developed for this purpose. Applications are demonstrated for the EA (Figure 7) and MTBE - I (Figures 8 and 9) systems.

For a system with nc components, critical equations are written, according to the following notation presented by de Medeiros (1990):




Formulae 42 - 44 were presented, in a different context, by Heidemann and Khalil (1980) for calculation of critical points in non-reactive systems using CEOS models. In problems without chemical reactions, the phase rule indicates nc – 1 degrees of freedom, generated by the existence of nc + 1 independent coordinates (T, P, x) and two constraints (equations 42 and 44, since equation 43 is used only for the calculation of vector u , with size nc x 1). Thus, normally, nc – 1 independent coordinates of vector x are specified and equations 42 and 44 are solved for (T, P).

For critical problems in reactive systems, the phase rule indicates nc nr – 1 degrees of freedom, due to the existence of nr additional chemical equilibrium constraints. In an analogy with the previous case, nc – nr – 1 compositional coordinates must be specified, and equations 42 and 44 and the chemical equilibrium constraints are then solved for (T, P) and the remaining compositional degrees of freedom. These specifications can be set using the nc – nr – 1 UD independent coordinates.

Algorithm 2 was implemented for determination of critical reactive coordinates (T, r) - RCP problem. This algorithm is based on the well-known non-reactive critical point calculation method proposed by Heidemann and Khalil (1980), which uses a double-loop structure of convergence for (), the phase composition given in molar fractions. In the present case, only the transformed phase composition () is specified (instead of molar fractions), whereas T and the molar densities vector () are sought numerically in order to satisfy equations 42, 43 and 44 and the chemical equilibrium constraints. Some characteristics of our algorithm are (i) an external iterative loop for the molar phase density () via Secant update; (ii) an internal structure for simultaneous resolution of chemical equilibrium equations, specification equations (in terms of the independent UD coordinates), the singularity condition for the Hessian matrix of Helmholtz free energy (equation 42) and the equation for temperature and molar densities vector ; (iii) compositional specification via independent UD coordinates () in order to match the nc - nr - 1 degrees of freedom; and (iv) compositional state described by molar densities.


Algorithm 2: Critical point calculation in reactive systems - RCP problem

1. are estimated and is specified;
2. l = 0 (counter);
3. are estimated;
4. System

is solved by the Newton-Raphson method for .

5. Image307.gif (1004 bytes) is solved for u;
6. is calculated;
7. is estimated, by the Secant method:

(If l = 0, is used).

8. Convergence is tested: if , stop.;
9. l = l + 1;
10. Steps 3-9 are repeated.


Figure 8 presents the bubble RVLE surfaces at 1 bar and 10 bar and the RCP surface (obtained using Algorithm 2) for the MTBE - I system using the vdW equation. Figure 9 presents the same objects predicted with the SRK equation. Algorithm 2 is also employed for the EA reactive system, whose RCP surface (using the SRK equation) is depicted in Figure 7.


Figure 8: Bubble RVLE surfaces and RCP surface for the MTBE -I system with the vdW equation.


Figure 9: Bubble RVLE surfaces and RCP surface for the MTBE -I system with the SRK equation.



In this work, phase rule calculations were performed for reactive systems under chemical equilibrium, using the UD coordinates (1995a). A specification formalism, suitable for conventional equations of state, is presented in order to generate reactive bubble, dew and azeotropic loci (i.e., first-order loci), and results were obtained for two example systems with industrial interest with two equations of state.

The generation of transformed molar enthalpy loci for saturated reactive systems revealed some evidence that the reactive McCabe-Thiele approximation (Platt, 1997) could be used in the design of reactive distillation operations with the MTBE - II system. On the other hand, a typical esterification (and more polar) system (EA) seems less likely to be described by this approximation.

With CEOS, the problem of determining reactive critical coordinates was also approached and RCP loci (i.e., third-order loci) were determined for the MTBE - I and EA systems. A fast and robust algorithm for determination of critical reactive states is presented.



Gustavo Mendes Platt is grateful to FAPERJ for its financial support of this research.



Latin letters

Aij binary interaction parameter of the Wilson model (K)

A, B, C, D constants for evaluation of ideal gas specific heat at constant pressure

AV Helmholtz free energy per unit of volume (J/m3)

â vector of activities

a CEOS attractive parameter

b CEOS covolume parameter

c cubic form of Helmholtz free energy per volume

specific heat at constant pressure for ideal gas (J/gmol.K)

CPL specific heat for liquids (J/gmol.K)

det determinant

E parameters for chemical equilibrium constant

CEOS cubic equation of state (vdW or SRK equation)

f fugacity (bar)

G Gibbs free energy (J/gmol)

k vapor-liquid equilibrium factor

kij interaction parameter of the SRK equation

H enthalpy (J/gmol)

enthalpy of formation as ideal gas at 298.15 K (J/gmol)

matrix defined in Table 9

M,m constant of reactive McCabe-Thiele approximation (J/gmol)

n constant of reactive McCabe-Thiele approximation (J/gmol)

nac number of additional constraints

nc number of components

np number of phases

nr number of equilibrium chemical reactions

P pressure (bar)

PvA, PvB, constants for vapor pressure (Wagner PvC, PvD equation)

pvA, pvB, constants for vapor pressure (Antoine pvC equation)

selection matrix

R universal gas constant (J/gmol.K)

selection matrix

T temperature (K)

u vector defined in equation 43

molar volume (m3/gmol)

X vector of UD coordinates for liquid phase

x vector of molar fractions for liquid phase

Y vector of UD coordinates for vapor phase

y vector of molar fractions for vapor phase


Greek letters

acentric factor

q dimensionless time

m chemical potential

D difference operator

f fugacity coefficient

gradient operator

Lij parameter of Wilson model

standard Gibbs free energy of reaction (J/gmol)

stoichiometric matrix

x vector defined in Table 9

lvap latent heat (J/gmol)

r molar density (gmol/m3)



a relative to independent UD  coordinates

C critical property

E excess property

G gas

i relative to component "i"

k relative to reference components (UD coordinates)

l, L liquid phase

R reduced property

React reactive

sat saturation condition

t transpose operator

v, V vapor phase



Aznar, M. and Telles, A. S., Cadernos de Termodinâmica (Tomo I). Departamento de Engenharia Química/Escola de Química/Universidade Federal do Rio de Janeiro, Brazil (in Portuguese) (1993).        [ Links ]

Barbosa, D. and Doherty, M. F., The Influence of Equilibrium Chemical Reactions on Vapor-Liquid Phase Diagrams, Chem. Eng. Sci., 43, p. 529 (1988a).        [ Links ]

Barbosa, D. and Doherty, M. F., The Simple Distillation of Homogeneous Reactive Mixtures, Chem. Eng. Sci., 43, p. 541 (1988b).        [ Links ]

Barbosa, D. and Doherty, M. F., Design and Minimum-Reflux Calculations for Single-feed Multicomponent Reactive Distillation, Chem. Eng. Sci., 43, p. 1523 (1988c).        [ Links ]

de Medeiros, J. L., Equilíbrio de Fases de Misturas Fluidas Próximo a Pontos Tricríticos. D. Sc. thesis (in Portuguese), PEQ/COPPE/ Universidade Federal do Rio de Janeiro, Brazil (1990).        [ Links ]

Graboski, M. S. and Daubert, T. E., A Modified Soave Equation of State for Phase Equilibrium Calculations. 1: Hydrocarbon Systems, Ind. Eng. Chem. Process Des. Dev., 17, p. 448 (1978).        [ Links ]

Heidemann, R. A. and Khalil, A. M., The Calculation of Critical Points, AIChE Journal, 26, p. 769 (1980).        [ Links ]

Ortega, J. M. and Rheinboldt, W. C., Iterative Solution of Nonlinear Equations in Several Variables. Academic Press (1970).        [ Links ]

Platt, G. M., Análise de Processos de Destilação Reativa. M.Sc. thesis (in Portuguese), TPQB/EQ/Universidade Federal do Rio de Janeiro, Brazil (1997).        [ Links ]

Reid, R. C., Prausnitz, J. M. and Poling, B. E., The Properties of Gases and Liquids. McGraw-Hill International Editions (1987).        [ Links ]

Smith, J. M., Van Ness, H. C. and Abbott, M. M., Introduction to Chemical Engineering Thermodynamics. McGraw-Hill International Editions (1996).        [ Links ]

Soave, G. S., Equilibrium Constants from a Modified Redlich-Kwong Equation of State, Chem. Eng. Sci., 27, p. 1197 (1972).        [ Links ]

Thiel, C., Sundmacher, K. and Hoffmann, U., Residue Curve Maps for Heterogeneously Catalysed Reactive Distillation of Fuel Ethers MTBE and TAME, Chem. Eng. Sci., 52, p. 1993 (1997).        [ Links ]

Ung, S. and Doherty, M. F., Vapor-liquid Equilibrium in Systems with Multiple Chemical Reactions, Chem. Eng. Sci., 50, p. 23 (1995a).        [ Links ]

Ung, S. and Doherty, M. F., Calculation of Residue Curve Maps for Mixtures with Multiple Chemical Reactions, Ind. Eng. Chem. Res., 34, p. 3195 (1995b).        [ Links ]

Walas, S. M., Phase Equilibria in Chemical Engineering. Butterworth Publishers (1985).        [ Links ]

Yu, J. and Eser, S., Determination of Critical Properties (Tc,Pc) of Some Jet Fuels, Ind. Eng. Chem. Res., 34, p. 404 (1995).        [ Links ]

Yu, J. and Eser, S., Thermal Decomposition of C10-C14 Normal Alkanes in Near-Critical and Supercritical Regions: Product Distributions and Reaction Mechanisms, Ind. Eng. Chem. Res., 36, p. 574 (1997a).        [ Links ]

Yu, J. and Eser, S., Kinetics of Supercritical-Phase Thermal Decomposition of C10-C14 Normal Alkanes and Their Mixtures, Ind. Eng. Chem. Res., 36, p. 585 (1997b).        [ Links ]

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License