Prediction of Tricritical Points in Ternary Mixtures by Global

Article pubs.acs.org/IECR

Prediction of Tricritical Points in Ternary Mixtures by Global Optimization Nélio Henderson,*,†,§ Raimundo A. Rodrigues, Jr.,†,§ and Wagner F. Sacco‡,§ †

Instituto Politécnico, Universidade do Estado do Rio de Janeiro, 28625-570, Nova Friburgo, Rio de Janeiro, Brazil Instituto de Engenharia e Geociências, Universidade Federal do Oeste do Pará, 68135-110, Santarém, Pará, Brazil § Thermodynamics and Optimization Group, Universidade do Estado do Rio de Janeiro, 28625-570, Nova Friburgo, Rio de Janeiro, Brazil ‡

ABSTRACT: Tricritical points of mixtures are characterized by thermodynamic coordinates where three coexisting phases become identical. Here, we formulate and solve for the first time the tricritical-point calculation as a global optimization problem. Using finite-difference formulas and the ordinary criticality criteria, we develop novel numerical approximations for tricriticality conditions associated with the derivatives of the fourth and fifth orders with respect to the compositions, obtained from the Gibbs tangent plane criterion for phase stability. Such approximations simplified the numerical implementations of these two conditions, avoiding propagation errors related to floating-point arithmetic without requiring a code with excessive numerical precision, which usually leads to high computational costs. This paper deals only with ternary mixtures. Thus, the present formulation is described as a minimization problem with strict inequality constraints, having four variables (two independent mole fractions, temperature, and pressure). Using the Peng−Robinson equation of state together with the classical one-fluid van der Waals mixing rule, we observe that the novel procedure is able to determine the tricritical coordinates in a set of ternary mixtures formed by alkanes or alkanes and light gases, which were previously studied by other authors and generally have complex multiphase behavior. The proposed method is compared to previous ones, and the results presented here compare favorably with the experimental values, which are available in the literature. In addition, we present an unprecedented prediction using the Peng−Robinson equation. We also show an example which leads to a no prediction of the tricritical point of a polar mixture, whose prediction (with this equation of state) was considered possible in a previous work.

1. INTRODUCTION Tricritical points of mixtures are characterized by thermodynamic coordinates where three coexisting phases become identical. In the 1960s, Radyshevskaya et al.1 and Krichevskii et al.2 were the first to publish reports of experimental realizations of tricritical behavior in mixtures. These authors observed the existence of tricritical points in mixtures with three and four components. A decade later, Windom3 and Griffiths and Windom4 observed a variety of mixtures which have a similar behavior. After that, Scott,5 Knobler and Scott,6 Windom and Sundar,7 and Lukes,8 among others, also analyzed mixtures subject to tricritical phenomena, including light gases and hydrocarbons studied by Kohn and co-workers.9−11 Mathematical formulations for the tricriticality criteria can be found in the works by Bartis,12 Mistura,13 and Michelsen and Heidemann,14 and in a recent article by Henderson et al.15 From the Gibbs phase rule, it is known that at least one component is necessary for the occurrence of an ordinary critical point, when two phases simultaneously become identical. Similarly, the Gibbs phase rule can be used to predict the minimum number of components in a mixture which has a tricritical point.16 In fact, according to this rule, an r-component mixture consisting of three phases in equilibrium has F = r − 3 + 2 degrees of freedom. Thus, when the three phases become identical, we must decrease F by 2, obtaining F = r − 3. Since F cannot be negative, then a tricritical point requires a minimum of r = 3 components, which occurs when F = 0. In this case, the tricritical point is said to be invariant, because it is the only © 2014 American Chemical Society

tricritical point of the referred ternary mixture. A mixture with three components can have a tricritical point and another different ordinary critical point, but such a mixture will have at most one single tricritical point. Furthermore, if, for a ternary mixture, there are two different points satisfying the tricriticality criteria, then more than three phases simultaneously become identical. As shown by Michelsen and Heidemann,14 for ternary mixtures the tricriticality criteria provide four nonlinear equations with four variables: two independent mole fractions, temperature, and pressure. Following Michelsen and Heidemann,14 the computation of tricritical points in ternary mixtures was specially treated by Kohse and Heidemann,17,18 using cubic equations of state. More recently, Arce and Aznar19 used the perturbed chain statistical associating fluid theory (PC-SAFT) equation of state20,21 to calculate tricritical points. All these authors employed numerical strategies based on the Newton method for nonlinear equations, whose convergence is tied to the location of the initial guess. Extending the methodology proposed by Henderson et al.22 for the prediction of ordinary critical points, here we calculate tricritical points of ternary mixtures via global optimization. This approach allows the use of robust global optimization Received: Revised: Accepted: Published: 4931

November 6, 2013 January 16, 2014 February 26, 2014 February 26, 2014 dx.doi.org/10.1021/ie403763a | Ind. Eng. Chem. Res. 2014, 53, 4931−4939

Industrial & Engineering Chemistry Research

Article r−1 r−1 r−1 r−1

algorithms that are easy to implement and show the ability to determine more than one solution, without depending on the localization of an initial guess and without calculation of derivatives of the objective function. This paper is organized as follows. In section 2 we summarize the tricriticality criteria obtained from Gibbs tangent plane criterion for phase stability and formulate the tricritical-point calculation as a global optimization problem. In the same section, we develop novel numerical approximations for tricriticality conditions associated with the derivatives of the fourth and fifth orders. This approach simplifies the numerical implementations of these two conditions, avoiding propagation errors related to floating-point arithmetic without requiring an excessive numerical precision, which usually leads to high computational costs. In section 3 we describe the optimization method used here, a swarm-based method originally developed for unconstrained global minimization, which was adapted to work with strict inequality constraints of a multimodal objective function. In section 4 we present the results of this work, and its comparisons with previous ones and with experimental data. For all mixtures considered here, we use the Peng−Robinson equation of state together with the classical one-fluid van der Waals mixing rule.23

C(T , P , z) =

i=1 j=1 k=1 l=1

D(T , P , z) =

(1)

i=1 j=1

r−1 r−1 r−1

B (T , P , z ) =

∑∑∑ i=1 j=1 k=1

∂ 3d(T , P , z) uiujuk ∂xi ∂xj ∂xk

A (T , P , z ) = 0

(6)

B (T , P , z ) = 0

(7)

C(T , P , z) = 0

(8)

D(T , P , z) = 0

(9)

∂d(T , P , x) = [μi (x) − μi (z)] − [μr (x) − μr (z)], ∂xi ∀ i = 1, ..., r − 1

(10)

Then, using eqs 1 and 10, and considering x = z, we have d(T,P,z) = 0 and ∇d(T,P,z) = 0. Consequently, given ε > 0, we can write the following Taylor’s expansion of d around z, from the term of order ε2:

A(T , P , z) = ∇2 d(T , P , z) ·u 2 ∂ 2d(T , P , z) uiuj ∂xi ∂xj

(5)

More specifically, u is an unitary eigenvector of the Hessian matrix of the tangent-plane distance function evaluated at the point (T,P,z), associated with the smallest eigenvalue of this matrix. In addition, the tricritical point (T,P,z) must be globally stable; i.e., given T, P, and z, there must occur d(T,P,x) ≥ 0, for all x ∈ Ω. To a systematic proof of these statements, see the recent paper by Henderson et al.15 2.2. Analytical Expressions and Approximations. In the context of this work, the prediction of tricritical points requires the calculation of the partial derivatives indicated in eqs 2−5, which can be described analytically or using numerical approximations. As mentioned before, in the present article the chemical potentials (and hence the tangent-plane distance function) are modeled from the Peng−Robinson equation of state with the classical one-fluid van der Waals mixing rule.23 Using this thermodynamic model, here the second order partial derivatives are computed analytically, as summarized by Henderson et al.22 Consequently, this approach will lead to the analytical calculation of the quadratic form described in eq 2. The forms in eqs 3−5 are estimated by finite-difference formulas. As shown by Henderson et al.,22 we note that the vector ∇d(x) = (∂d(x)/∂xi), the gradient of d, is determined analytically by

where T and P are, respectively, the temperature and the pressure of the mixture, and μj represents the chemical potential of component j = 1, ..., r, which also depends on T and P, i.e., μj(x) = μj(T, P, x). According to the so-called Gibbs tangent plane (GTP) criterion for phase stability,24 the thermodynamic mixture with composition z is stable, at temperature T and pressure P, if and only if d(T,P,x) ≥ 0, for all x ∈ Ω. Tricriticality criteria for an r mixture with composition z ∈ Ω can be obtained by expanding the tangent-plane distance function in Taylor series around z.15 For example, in an expansion of the form d(T, P, z + uε) = (A/2!)ε2 + (B/3!)ε3 + (C/4!)ε4 + (D/5!)ε5 + ... arise partial derivatives of higher order of the function d. Here, to cope with these derivatives, we consider the following representations:

∑∑

∂ 5d(T , P , z) ∂xi ∂xj ∂xk ∂xl ∂xm

where the vector u = (u1, ..., un−1) ∈ r − 1, a function of T, P, and z, will be defined as follows. If (T,P,z), with z ∈ Ω, is a tricritical point of an r-component mixture, then, from the GTP criterion, it can be shown15 that there is a vector u = u(T,P,z) ∈ r − 1, with u·u = 1, such that

i=1

=

(4)

uiujukulum

∑ xi{[μi (x) − μi (z)] − [μr (x) − μr (z)]}

r−1 r−1

∑∑∑∑ ∑ i=1 j=1 k=1 l=1 m=1

r−1

+ [μr (x) − μr (z)]

∂ 4d(T , P , z) uiujukul ∂xi ∂xj ∂xk ∂xl

r−1 r−1 r−1 r−1 r−1

2. THEORY 2.1. Tricriticality Conditions from the Gibbs Tangent Plane Criterion. Let Ω = {x = (x1,...,xr−1) ∈ r − 1; 0 < xi < 1 and ∑r−1 i=1 xi < 1} be the compositional space of r-component mixtures, where the mole fraction xr is obtained from xr = 1 − ∑r−1 i=1 xi. Given a mixture with composition z = (z1,..., zr−1) ∈ Ω, the called tangent-plane distance function24 can be written as22 d (T , P , x ) =

∑∑∑∑

(2)

d(T , P , z − εu) =

A 2 B C D ε − ε 3 + ε 4 − ε 5 + ... 2! 3! 4! 5! (11)

Performing the differentiation of eq 11 with respect to the parameter ε, we obtain

(3) 4932

dx.doi.org/10.1021/ie403763a | Ind. Eng. Chem. Res. 2014, 53, 4931−4939


Article

d B C D d(T , P , z − εu) = Aε − ε 2 + ε 3 − ε 4 + ... dε 2! 3! 4!

D=

(12)

On the other hand, using the chain rule, we have d d(T , P , z − εu) = −∇d(T , P , z − εu) ·u dε

where the partial derivatives of second order ∂2d(T,P,z+εu)/∂xi ∂xj and ∂2d(T,P,z−εu)/∂xi ∂xj are also calculated analytically.22 2.3. Some Simplifications and the Optimization Problem. In this section, we describe the optimization problem formulated for the calculation of tricritical points of ternary mixtures. But before that, combining the approximations in eqs 17 and 23 with the ordinary criticality criteria indicated in eqs 6 and 7, we develop a new approach that simplifies the implementation of the tricriticality conditions shown in eqs 8 and 9. Using ε of the order of magnitude of 10−6, in practice we had the opportunity to observe that, during iterative processes, the implementation of the formulas 17 and 23 (whose denominators reached the order of 10−18) were affected by numerical problems associated with the floating-point arithmetic, even when encoded in double precision, causing (with relative frequency) the divergence of the sequence generated by the numerical method. This problem is not new. In fact, in order to avoid such difficulties, Arce and Aznar,19 for example, used a FORTRAN program with a numerical precision of 32 digits, in the implementation of the finite-difference formulas relating to higher-order partial derivatives of molar Gibbs free energy of mixtures, which were used as mathematical criteria for tricritical points. Here, we do not use the excessive numerical precision of 32 digits, because such an approach is computationally expensive. To overcome this numerical problem, initially, from eq 17, note that we can consider the following approximation for the tricriticality condition in eq 8.

(13)

Combining eqs 12 and 13, we arrive at the relation ∇d(T , P , z − εu) ·u = −Aε +

B 2 C D ε − ε 3 + ε 4 − ... 2! 3! 4! (14)

Using an entirely analogous procedure, we can also write ∇d(T , P , z + εu) ·u = Aε +

B 2 C D ε + ε 3 + ε 4 + ... 2! 3! 4! (15)

By adding eq 15 to eq 14, we get the following second-order approximation for the form in eq 3: B=

1 [∇d(z + εu) + ∇d(z − εu)]·u + O(ε 2) ε2

(16)

Now, to approximate the form shown in eq 4, we subtract eq 14 from eq 15, obtaining a formula with a leading error of second order, given by C=

3 6 [∇d(z + εu) − ∇d(z − εu)]·u − 2 A + O(ε 2) 3 ε ε (17)

Note that the vectors ∇d(z + εu) and ∇d(z − εu) shown in eqs 16 and 17 are computed by the analytical formulas described in eq 10.22 To determine an approximation for the form shown in eq 5, initially we perform the differentiation of eq 12, obtaining

1 [∇d(z + εu) − ∇d(z − εu)]·u = 2A ε

d2 C D d(T , P , z − εu) = A − Bε + ε 2 − ε 3 + ... 2 2! 3! dε

1 2 [∇ d(T , P , z + εu) ·u 2 − ∇2 d(T , P , z − εu) ·u 2] = 2B ε

Again, using the chain rule, one can show that

(25)

Since at a tricritical point the ordinary conditions A = B = 0 also occur, then, replacing these zero values in eqs 24 and 25, we obtain C̃ (T,P,z) = D̃ (T,P,z) = 0, where the functions C̃ and D̃ are defined, respectively, by the left sides of eqs 24 and 25. In other words, in the present work the mathematical criteria for a tricritical point will be approximated by the simplified relations:

(19)

where r−1 r−1

∇2 d(T , P , z − εu) ·u 2 =

∑∑ i=1 j=1

∂ 2d(T , P , z − εu) uiuj ∂xi ∂xj (20)

A(T , P , z) = ∇2 d(T , P , z) ·u 2 = 0

Thus, from eqs 18 and 19, we have

1 [∇d(z + εu) + ∇d(z − εu)]·u ε2 =0 1 C̃(T , P , z) = [∇d(z + εu) − ∇d(z − εu)]·u ε =0 1 D̃ (T , P , z) = [∇2 d(T , P , z + εu) ·u 2 ε − ∇2 d(T , P , z − εu) ·u 2] =0 B (T , P , z ) =

C D ∇ d(T , P , z − εu) ·u = A − Bε + ε 2 − ε 3 + ... 2! 3! 2

2

(21)

Similarly, it can be shown that ∇2 d(T , P , z + εu) ·u 2 = A + Bε +

(24)

Similarly, using the formula shown in eq 23, the tricriticality condition summarized in eq 9 can be approximated by

(18)

d2 d(T , P , z − εu) = ∇2 d(T , P , z − εu) ·u 2 dε 2

3 2 [∇ d(T , P , z + εu) ·u 2 − ∇2 d(T , P , z − εu) ·u 2] ε3 6 − 2 B + O(ε 2) (23) ε

C 2 D ε + ε 3 + ... 2! 3! (22)

Finally, subtracting eq 21 from eq 22, we obtain the following second-order approximation by finite difference: 4933

(26)



Article

Now, the finite-difference formulas shown in eq 26 are easily implemented using double precision, without causing propagation errors related to floating-point arithmetic. For ternary mixtures (r = 3), eq 26 describes a nonlinear system with four equations and four variables, T, P, and the compositions z1 and z2, which constitute the vector z = (z1, z2) ∈ Ω. Thus, note that the mole fraction z3 is obtained from z3 = 1 − (z1 + z2). Let (Tmin, Tmax) and (Pmin, Pmax) be open intervals where we can find, respectively, the temperature and the pressure of a tricritical point with composition z = (z1,z2). Extending the previous work by Henderson and co-workers22 for the calculation of ordinary critical points, here we formulate the tricritical-point calculation in ternary mixtures as the following global optimization problem:

population-based algorithm, where the population is called a swarm and each individual is called a particle. Here, the PSO algorithm is addressed to global optimization problems having the following form, which summarizes the minimization problem formulated in eq 27: min f = f (x1 , x 2 , x3 , x4)

subject to a1 < x1 < b1

a 2 < x 2 < b2 a3 < x3 < b3 a4 < x4 < b4 − x3

where aj and bj are, respectively, lower and upper bounds on the variables. PSO uses a swarm of np particles, where np is the population size. Here, we consider np = 100. At a given instant of time t, each particle is associated with one vector position, denoted by xi(t) = (xi1(t), ..., xi4(t)), for all i = 1, ..., np. In the numerical context, PSO is an iterative method, where t is the current iteration (or cycle). Thus the new position of the ith (i = 1, ..., np) particle in the iteration t + 1 is given by

min f (T , P , z) = A(T , P , z)2 + B(T , P , z)2 + D̃ (T , P , z)2 + C̃(T , P , z)2

(28)

(27)

subject to

Tmin < T < Tmax Pmin < P < Pmax z1min < z1 < z1max

xji(t + 1) = xji(t ) + vji(t + 1);

z 2min < z 2 < 1 − z1

j = 1, ..., 4

(29)

Here, vi(t+1) = (vi1(t+1), ..., vi4(t+1)) is a vector in 4 that represents the velocity of the particle xi at time t + 1. To control possible divergences (“explosions”) of the swarm, Clerc and Kennedy26 proposed a robust version of the PSO method, called PSO with constriction factor (PSOCF). In this version, the velocity vector vi(t+1) is updated according to the equation

where z1min and z1max are, respectively, the lower and upper bounds on the variable z1, and z2min is the lower bound on the variable z2. The inequality z2 < 1 − z1 ensures that z3, calculated by z3 = 1 − (z1 + z2), is positive. A global minimizer (T,P,z) of this bound-constrained optimization problem, where the objective function f takes the value 0, is a solution of the nonlinear system shown in eq 26. In addition, such a global minimizer will be considered a tricritical point of the ternary system if, for this composition z, d(T,P,x) ≥ 0, for all x ∈ Ω. Here, the finite difference approximations considered above are computed with ε = 10−6∥z∥, where z is the current composition, and the eigenvector calculations of the Hessian matrix of d are performed using the Jacobi method.25

vji(t + 1) = χ[vji(t ) + ϕ1r1j(yji (t ) − xji(t )) + ϕ2r2j(xbest j(t ) − xji(t ))];

j = 1, ..., 4

(30)

where ϕ1 and ϕ2 are positive constants, and r1j and r2j are random numbers uniformly distributed in [0, 1]. yij(t) is the jth component of the vector yi(t), the best previous position encountered by the ith particle, and xbestj(t) is the jth component of the vector xbest(t), the best previous position among all the particles of the swarm. The parameter χ is the constriction factor. Using a theoretical analysis, Clerc and Kennedy26 were able to demonstrate that this constriction factor must satisfy the equation

3. THE NUMERICAL METHOD During the development of this work, we could observe that the calculation of tricritical points via optimization is a difficult problem that requires robust methods for global minimization of a multimodal objective function, restricted to an open set with lower and upper bounds and having one additional strict inequality constraint. In this paper, we use a modification of a robust version of the particle swarm optimization (PSO) method proposed by Clerc and Kennedy.26 As originally proposed by Eberhart and Kennedy,27 the PSO method is based on the collaborative paradigm and social learning metaphor, applied to unconstrained optimization problems. Initially inspired by the dynamics of fish schools and bird flocks, it balances the learning and experience of a particle and its ability in taking advantage of the experience of the swarm in a social context, where each individual learns from its own experience, as from the experience of the group, balancing exploration and exploitation. This metaheuristic is a

2

χ= |2 − ϕ −

ϕ2 − 4ϕ |

(31)

where ϕ ≡ ϕ1 + ϕ2 > 4.0

(32)

Typically, ϕ is set as 4.1 and the constriction factor χ is then 0.729. Here, we consider ϕ1 = ϕ2 = 2.05. PSOCF is a stochastic method for unconstrained global optimization. In order to adapt this method to the constrained problem shown in eq 28, in this work we modify PSOCF, forcing, in a suitable manner, the sequence generated by the PSOCF method remain within the feasible set. This is done as follows. For all j = 1, 2, 3, after the iteration that is described by 4934



Article

eqs 29 and 30, if xij(t+1) ≥ bj, then we reset the velocity vij(t+1), making vji(t + 1) = bj − xji(t );

Similarly, if

xij(t+1)

j = 1, 2, 3

(33)

≤ aj, then we redefine the velocity as

vji(t + 1) = aj − xji(t );

j = 1, 2, 3

(34)

After such changes, we consider xji(t + 1) = xji(t ) + λijvji(t + 1);

j = 1, 2, 3

(35)

where λij are random numbers uniformly distributed in the interval (0,1). To impose the inequalities a4 < xi4 < b4 − xi3, we write these constraints as a4 < x4i < β i

(36)

where the number β = b4 − is obtained after the calculation of xi3. Therefore, after using the updates shown in eqs 29 and 30, if xi4(t+1) ≥ βi, or if xi4(t+1) ≤ a4, then we reset the velocity vi4(t+1) as vi4(t+1) = βi − xi4(t), or as vi4(t+1) = a4 − xi4(t), and then we consider xi4(t+1) = xi4(t) + λi4vi4(t+1), where, for all i, λi4 are random numbers uniformly distributed in the interval (0,1). The steps of the PSOCF with the modifications mentioned above (here called PSOCF/M) are shown in algorithm 1 in Figure 1. In algorithm 1 (Figure 1), the iterative process terminates only when the maximum number of cycles is reached. Here, we use 10 000 cycles. Our computational experiments also indicate that the tricritical-point calculation with algorithm 1 has considerable sensitivity to the location of the initial particles. To overcome this problem, in the present work the initial population of this algorithm is generated by the Sobol sequences,28 which are lowdiscrepancy sequences.29 This is done by considering {ζi1, ζi2, ζi3, ζi4)}i=1,2,... as a Sobol sequence in 1 , where 1 = (0, 1)4 is the open unitary cube in 4 . In the present work, we use a code for Sobol sequences available in the Numerical Recipes collection.30 i

xi3

4. RESULTS In order to test and compare the proposed methodology, we consider five ternary mixtures previously analyzed by Kohse and Heidemann. 17,18 Among all the ternary mixtures considered by these authors, we study only those that have the following properties: (i) the mixture can be appropriately modeled by the Peng−Robinson (PR) equation of state together with the classical one-fluid van der Waals mixing rule, (ii) the tricritical coordinates have been previously obtained experimentally, and (iii) the tricritical temperature and tricritical pressure were also calculated by Arce and Aznar,19 using the PC-SAFT (perturbed chain statistical associating fluid theory) equation of state. In addition, we consider the study of two other mixtures. In fact, we present an unprecedented prediction using the Peng− Robinson equation, for the ternary mixture comprising methane, n-pentane, and n-octane. We also show an example which leads to a no prediction of the tricritical point of a polar mixture, whose prediction (with this equation of state) was considered possible by Kohse and Heidemann.18

Figure 1. Algorithm 1: steps of the PSOCF with the modifications as cited in the text (PSOCF/M).

The input data for each pure component required by the PR equation, namely the critical temperatures, critical pressures, and acentric factors, were taken from Kohse’s thesis.17 For each tested ternary mixture, in the mixing rule we use the nonzero binary interaction parameters also indicated by Kohse.17 Our numerical experiments have shown that the objective function of the optimization problem described in eq 27 usually has different minimizers. Thus, to obtain the tricritical point, we use an appropriate selection. First, for each problem, numerical experiments are repeated 100 times, using distinct random seeds for the random number generating developed by Matsumoto and Nishimura.31 Second, using the GTP criterion for phase stability, we test the stability of each calculated minimizer. The unstable points are discarded. Finally, among the remaining minimizers, the selected tricritical point is the one with smaller value of the objective function. To perform 4935



Article

method are consistent with the experimental data.9 Although the calculated composition z3 shows a large error in comparison with the experimental value, we note that the other tricritical coordinates calculated by Kohse and Heidemann17,18 are in agreement with these experimental data. The same happens with the tricritical temperature and tricritical pressure calculated from the PC-SAFT equation.19 In this example, the objective function used here reached the order of 10−7. 4.3. Mixture 3: CH4 (1)/n-C4H10 (2)/n-C8H18 (3). This mixture was also analyzed experimentally by Kohn and Luks.9 In this example we consider Tmin = 200 K, Tmax = 300 K, Pmin = 6000 kPa, Pmax = 7000 kPa, and z1max = 0.95. As can be noted from Table 3, for this mixture the tricritical point calculated by the PSOCF/M method is in agreement with

stability tests, this selection requires a reliable method. We use the PJSA (projected simulated annealing) algorithm developed by Henderson et al.32 Here, in all examples, we use z1min = z2min = 5 × 10−4. 4.1. Mixture 1: CH4 (1)/C2H6 (2)/n-C8H18 (3). This mixture was experimentally studied by Kohn and Luks.9 In this example, we consider Tmin = 200 K, Tmax = 300 K, Pmin = 6500 kPa, Pmax = 7500 kPa, and z1max = 0.85. Table 1 shows the experimental tricritical coordinates, the results obtained in this work and those calculated by Kohse and Table 1. Comparison of Calculated and Experimental Tricritical Points for the Mixture CH4/C2H6/n-C8H18

exptl9 this worka Kohse and Heidemann18 a Arce and Aznar19 b a

T (K)

P (MPa)

z2

z3

f

222.0 222.0 234.0

6.991 6.774 7.556

0.200 0.197 0.281

0.035 0.006 0.014

− 2.5 × 10−9 −

223.4

6.924

−

−

−

Table 3. Comparison of Calculated and Experimental Tricritical Points for the Mixture CH4/n-C4H10/n-C8H18

9

b

exptl this worka Kohse and Heidemann18 a Arce and Aznar19 b

PR equation. PC-SAFT equation.

T (K)

P (MPa)

z2

z3

f

207.3 205.9 203.2

6.520 6.624 6.047

0.0620 0.0829 0.0486

0.00960 0.00230 0.00069

− 1.3 × 10−7 −

206.8

6.472

−

−

−

Heidemann,17,18 both using the PR equation with nonzero binary interaction parameters, and the tricritical temperature and tricritical pressure calculated by Arce and Aznar,19 from the PC-SAFT equation. Unfortunately, this work of Arce and Aznar does not provide the calculated tricritical compositions. From Table 1, we can note that the tricritical temperature calculated by the present method coincides with the experimental tricritical temperature, and the tricritical pressure is relatively close to the experimental tricritical pressure. On the other hand, we can note that the tricritical temperature and tricritical pressure calculated by the method proposed by Kohse and Heidemann17,18 do not have the same degree of accuracy with regard to the experimental data.9 Furthermore, our results, which were obtained with the PR equation together with the classical mixing rule, are competitive with those calculated by Arce and Aznar,19 using the modern PC-SAFT equation. However, we can note that the calculated tricritical composition z3 = 0.006 does not exhibit this same accuracy. The objective function used here reached the order of 10−9, which may be considered significant, taking into account the simplifications and approximations shown in eq 26. 4.2. Mixture 2: CH4 (1)/C3H6 (2)/n-C8H18 (3). This mixture was also analyzed experimentally by Kohn and Luks.9 Here, the values of Tmin, Tmax, Pmin, Pmax, and z1max are the same as in the previous example. With exception of the calculated tricritical composition z3, Table 2 shows that the results obtained using the present

the experimental data, including the value of the calculated tricritical composition z3. The same accuracy cannot be observed at the tricritical point calculated by Kohse and Heidemann,17,18 particularly with respect to the composition z3. On the other hand, the tricritical temperature and tricritical pressure that were calculated from the PC-SAFT equation19 have similar accuracies. Thus, these results allow us to conclude that (for this mixture) the loss of accuracy associated with the tricritical point obtained by the method of Kohse and Heidemann18 is not due to the use of the PR equation, together with the classical one-fluid van der Waals mixing rule. In fact, for this mixture of alkanes, the PR equation and the PCSAFT equation are both appropriate models for the calculation of tricritical points. In this example, the objective function f reached the order of 10−7. 4.4. Mixture 4: CH4 (1)/CO2 (2)/n-C8H18 (3). This mixture was also analyzed experimentally by Kohn and Luks.9 In this example we consider Tmin = 200 K, Tmax = 300 K, Pmin = 6500 kPa, Pmax = 7500 kPa, and z1max = 0.85. From Table 4, we can note that the tricritical point calculated by the PSOCF/M method is entirely in good agreement with the experimental data,9 and the tricritical temperature and

Table 2. Comparison of Calculated and Experimental Tricritical Points for the Mixture CH4/C3H6/n-C8H18

Table 4. Comparison of Calculated and Experimental Tricritical Points for the Mixture CH4/CO2/n-C8H18

exptl9 this worka Kohse and Heidemann18 a Arce and Aznar19 b a

T (K)

P (MPa)

z2

z3

f

215.3 215.2 215.1

7.113 7.188 7.122

0.106 0.141 0.115

0.1550 0.0119 0.0061

− 1.1 × 10−7 −

215.2

7.103

−

−

−

a

PR equation. bPC-SAFT equation.

exptl9 this worka Kohse and Heidemann18 a Arce and Aznar19 b

PR equation. bPC-SAFT equation.

a

4936

T (K)

P (MPa)

z2

z3

f

222.0 220.7 215.1

7.306 7.463 6.488

0.2990 0.2854 0.2270

0.0365 0.0144 0.0049

− 4.9 × 10−12 −

221.5

7.291

−

−

−

PR equation. bPC-SAFT equation. dx.doi.org/10.1021/ie403763a | Ind. Eng. Chem. Res. 2014, 53, 4931−4939


Article

tricritical pressure are also competitive with the corresponding values calculated by the PC-SAFT equation.19 Despite having been determined by the same thermodynamic model, the tricritical point calculated by Kohse and Heidemann17,18 does not have the same concordance. Thus, for this mixture of alkanes and CO2, this fact indicates that the differences shown previously by Arce and Aznar19 are not associated with the PR equation. In fact, probably this loss of accuracy is a problem intrinsic to the local search method used by Kohse and Heidemann.17,18 At the tricritical point calculated by this methodology, the objective function f reached the significant order of 10−12. 4.5. Mixture 5: CH4 (1)/N2 (2)/n-C4H10 (3). This mixture was experimentally studied by Kohn et al.11 In this example we consider Tmin = 170 K, Tmax = 200 K, Pmin = 6000 kPa, Pmax = 7000 kPa, and z1max = 0.85. Table 5 shows that the tricritical coordinates calculated using the PSOCF/M method compare favorably with the exper-

Table 6. Experimental and Calculated Tricritical Points for the Mixture CH4/n-C5H12/n-C8H18 exptl10 this worka a

exptl11 this worka Kohse and Heidemann18 a this worka a

P (MPa)

z2

z3

f

175.5 176.6 194.0

6.219 6.001 6.239

0.2350 0.3000 0.0965

0.1230 0.0120 0.0401

− 8.9 × 10−8 −

173.8

6.225

−

−

−

P (MPa)

z2

z3

f

6.10 5.60

0.060 0.095

0.0050 0.0014

− 2.9 × 10−7

PR equation.

Robinson cubic equation of state, it is possible to make a very reasonable numerical prediction for a tricritical point of this important ternary mixture, commonly found in the oil industry. These results were generated using Tmin = 150 K, Tmax = 400 K, Pmin = 5000 kPa, Pmax = 8000 kPa, and z1max = 0.95. 4.7. An Example Which Leads to a No Prediction. In the context of the present methodology, it is considered that a tricritical point cannot be estimated (or there is no such point) if the objective function described in eq 27 has not reached a sufficiently small value and/or if the calculated minimizers are unstable; i.e., such minimizers do not satisfy the Gibbs tangent plane criterion for phase stability. In this section, we consider an example where this occurs, more precisely, where all the calculated minimizers are always unstable. In one of the pioneering works concerning the study of tricritical phenomena, using laboratory experiments, Shvarts and Efremova33 determined one tricritical point for the ternary polar mixture comprising carbon dioxide, ethanol, and water. Employing the Peng−Robinson cubic equation of state, together with the classical one-fluid van der Waals mixing rule, Kohse and Heidemann18 calculated a point for this polar mixture, which exhibits a clear discrepancy concerning the experimental tricritical pressure value of Shvarts and Efremova;33 see Table 7. As shown in Table 7, using the PJSA algorithm developed by Henderson et al.,32 we were able to verify that this point calculated by Kohse and Heidemann18 is (globally) unstable. Thus, the point determined by these authors18 cannot be regarded as a tricritical point of the polar mixture CO2/ C2H5OH/H2O, predicted by the PR equation. This is explained by the characteristic of the stability test used by Kohse and Heidemann,18 which is just a local stability test not appropriate for analyzing tricritical points. Employing Tmin = 300 K, Tmax = 400 K, Pmin = 9000 kPa, Pmax = 12 000 kPa, and z1max = 0.95, after extensive testing with the proposed method, we verified that all the calculated minimizers are always unstable. On the basis of these results, we can observe that the Peng− Robinson cubic equation of state (equipped with the classical one-fluid van der Waals mixing rule) is unable to predict tricritical phenomena in the ternary polar mixture comprising carbon dioxide, ethanol, and water.

Table 5. Comparison of Calculated and Experimental Tricritical Points for the Mixture CH4/N2/n-C4H10 T (K)

T (K) 203.30 192.65

PR equation.

imental values of Kohn et al.11 The same can be said with respect to the tricritical temperature and tricritical pressure calculated from the PC-SAFT equation.19 However, the results generated by the method proposed by Kohse and Heidemann17,18 are flawed. In fact, the calculated tricritical temperature and all the calculated tricritical compositions show errors in comparison with the experimental data. Again, in view of the results obtained by the present methodology, we can say that, for this mixture of alkanes and N2, such failures are not a consequence of possible nonphysical predictions by the Peng− Robinson equation of state. At the tricritical point calculated by PSOCF/M method, the objective function f reached the order of 10−8. 4.6. An Unprecedented Prediction Using the PR Equation. When referring to the cubic equations of state, Kohse and Heidemann18 stated, “The inability of the model to locate a tricritical point for the methane + nitrogen + n-pentane system is due to the equations of state incorrectly predicting liquid−liquid immiscibility for the methane + n-pentane binary pair. This immiscibility only occurs with alkanes of carbon number of 6 or higher. This is also the reason that the models are unable to locate tricritical points for any systems in the series CH4 + nC5H12 + nCNH2N+2” (see p 1252 of ref 18). In this section, we show that this statement made by Khose and Heidemann18 is not an absolute truth. For this, using the PR equation of state, we investigate the existence of a tricritical point associated with a system just mentioned by these authors, the ternary mixture comprising CH4/n-C5H12/n-C8H18. Thus, employing the methodology developed here, we obtained the results shown in Table 6, which also displays the experimental data generated by Kohn and co-workers.10 As we can see, contrary to the above statement, using the Peng−

5. CONCLUSIONS Using tricriticality conditions that were previously deduced by these authors from a systematic approach,15 in the present work we formulate and solve the tricritical point calculation problem of thermodynamic mixtures as an optimization problem. In addition, using finite-difference formulas and the ordinary criticality conditions, we developed novel numerical approximations for tricriticality criteria associated with the forms of the fourth and fifth degrees. Such approximations simplified the numerical implementations of these two conditions, avoiding 4937



Article

Table 7. Experimental Tricritical Point for the Mixture CO2/C2H5OH/H2O and the Impossibility of Numerical Predictions Using the PR Equation exptl33 Kohse and Heidemann18 a this worka a

T (K)

P (MPa)

z2

z3

state

320.6 334.1 NTCPb

9.292 11.283 NTCPb

0.1500 0.0908 NTCPb

0.0200 0.0640 NTCPb

− unstable unstable

PR equation. bNTCP, no tricritical points.

■

propagation errors related to floating-point arithmetic without requiring high computational costs. This work was addressed to ternary mixtures. Thus, the optimization problem was introduced as a global minimization problem with four variables (two independent mole fractions, temperature, and pressure), having strict inequality constraints. Due to the nature of the tricritical-point calculation in ternary mixtures, we might note that the problem formulated here is a difficult minimization problem with linear constraints, which usually has several minimizers. In order to overcome such difficulties, and to avoid possible divergences of the iterative process, we employ here a robust version of the particle swarm optimization (PSO) method proposed by Clerc and Kennedy,26 which was modified to work with inequality constraints. To test this new methodology, five systems were initially considered, essentially mixtures of alkanes or alkanes and light gases, which can be appropriately modeled by the Peng− Robinson equation of state (together with the classical onefluid van der Waals mixing rule) used in our formulation. The tricritical coordinates of these systems were experimentally obtained by Kohn and co-workers,9,11 and previously calculated by Kohse and Heidemann,17,18 from the Peng−Robinson equation of state. Using the PC-SAFT equation of state, Arce and Aznar19 also calculated the tricritical temperature and tricritical pressure associated with such systems. In general, for these mixtures, the results presented here compare favorably with the experimental values,9,11 are superior to those obtained by Kohse and Heidemann,17,18 and are competitive with the tricritical temperatures and tricritical pressures calculated from the modern PC-SAFT equation.19 Finally, in addition, we showed an unprecedented prediction for the ternary mixture comprising methane, n-pentane, and noctane, and the no prediction of the tricritical point of the polar mixture CO2/C2H5OH/H2O, whose prediction was considered to be possible by Kohse and Heidemann.18 Thus, these facts allow us also to conclude that same failures observed in the previous results of Kohse and Heidemann17,18 are not (necessarily) consequences of possible nonphysical predictions by the Peng−Robinson equation of state. Probably, such failures are consequence of intrinsic limitations concerning the local search method used by these authors. The extension of the present approach to mixtures containing more than three components will be the subject of future works.

■

ACKNOWLEDGMENTS N.H. and W.F.S. gratefully acknowledge the financial support provided by CNPq (Conselho Nacional de Desenvolvimento ́ Cientifico e Tecnológico, Ministry of Science & Technology, Brazil). R.A.R., Jr., was supported by CAPES (Coordenaçaõ de ́ Superior, Brazil). The Aperfeiçoamento de Pessoal de Nivel research by N.H. has been carried out in the framework of Project PROCIENCIA-UERJ financed by FAPERJ.

■

REFERENCES

(1) Radyshevskaya, G. S.; Nikurashina, N. I.; Mertslin, R. V. Effect of Temperature on Equilibria or Three Liquid Phases in Fourcomponent Systems. J. Gen. Chem. (USSR) 1962, 32, 673−676. (2) Krichevskii, I. R.; Efremova, G. D.; Pryanikova, R. O.; Serebryakova, A. V. A Possible Case of Critical Phenomena. Russ. J. Phys. Chem. 1963, 37, 1046−1047. (3) Windom, B. Tricritical Points in Three- and Four-component Fluid Mixtures. J. Phys. Chem. 1973, 77, 2196−2200. (4) Griffiths, R. B.; Windom, B. Multicomponent-Fluid Tricritical Points. Phys. Rev. A 1973, 8, 2173−2175. (5) Scott, R. L. Multicritical Phenomena in Fluid Mixtures. In Proceedings of the Eighth Symposium on Thermophysical Properties, Vol. 1: Thermophysical Properties of Fluids; Sengers, J. V., Ed.; American Society of Mechanical Engineers: New York, 1982; pp 397−404. (6) Knobler, C. M.; Scott, R. L. Multicritical Points in Fluid Mixtures: Experimental Studies. Phase Transitions and Critical Phenomena; Domb, C., Lebowitz, J. L., Eds.; Academic Press: New York, 1984; Vol. 9, Chapter 2. (7) Windom, B.; Sundar, G. Some Theoretical Aspects of Critical Phenomena in Fluids. Fluid Phase Equilib. 1986, 30, 1−14. (8) Luks, K. D. The Occurrence and Measurement of Multiphase Equilibrium Behavior. Fluid Phase Equilib. 1986, 29, 209−224. (9) Kohn, J. P.; Luks, K. D. Liquid-Liquid-Vapor Equilibria in Cryogenic LNG Mixtures; Research Report RR-49; Gas Processors Association: Tulsa, OK, 1981. (10) Kohn, J. P.; Merrill, R. C.; Luks, K. D. Liquid-Liquid-Vapor Equilibria in Cryogenic LNG Mixtures. Phase II; Research Report RR60; Gas Processors Association: Tulsa, OK, 1982. (11) Kohn, J. P.; Merrill, R. C.; Luks, K. D. Liquid-Liquid-Vapor Equilibria in Cryogenic LNG Mixtures. Phase III-Nitrogen Rich Systems; Research Report RR-67; Gas Processors Association: Tulsa, OK, 1983. (12) Bartis, J. T. Thermodynamic Equations for Tri- or Third Order Critical Points. J. Chem. Phys. 1973, 59, 5423−5430. (13) Mistura, L. Critical Phases in Multicomponent Fluid Mixtures. J. Phys. A 1976, 9, 2139−2148. (14) Michelsen, M.; Heidemann, R. A. Calculation of Tri-Critical Points. Fluid Phase Equilib. 1988, 39, 53−74. (15) Henderson, N.; Sacco, W. F.; Rodrigues, R. A., Jr. A Deduction of the Multicriticality Conditions of Mixtures from the Gibbs Tangent Plane Criterion. Fluid Phase Equilib. 2013, 360, 305−308. (16) Zernike, L. General Considerations Concerning the Number of Virtual Phases. Recl. Trav. Chim. Pays-Bas 1949, 68, 585−594. (17) Kohse, B. F. The Calculation of Tricritical Points. M.S. Thesis, The University of Calgary, Canada, 1989. (18) Kohse, B. F.; Heidemann, R. A. Computation of Tricritical Points in Ternary Systems. AIChE J. 1993, 39, 1242−1256. (19) Arce, P.; Aznar, M. Computation and Modeling of Tricritical Phenomena in Ternary and Quaternary Mixtures Using the Perturbed

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +55 (22) 2533-2263. Fax: +55 (22) 2533-2263. Notes

The authors declare no competing financial interest. 4938



Article

Chain-Statistical Associating Fluid Theory. J. Supercrit. Fluids 2009, 49, 135−142. (20) Gross, J.; Sadowski, G. Application of Perturbation Theory to a Hard-Chain Reference Fluid: An Equation of State for Square-Well Chains. Fluid Phase Equilib. 2000, 168, 183−199. (21) Gross, J.; Sadowski, G. Perturbed-Chain, SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (22) Henderson, N.; Freitas, L.; Platt, G. M. Prediction of Critical Points: A New Methodology Using Global Optimization. AIChE J. 2004, 50, 1300−1314. (23) Peng, D.; Robinson, B. D. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (24) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. SPEJ, Soc. Pet. Eng. J. 1982, 22, 731−742. (25) Sewell, G. Computational Methods of Linear Algebra; Ellis Horwood: London, 1990. (26) Clerc, M.; Kennedy, J. The Particle Swarm-Explosion, Stability and Convergence in a Multidimensional Complexe Space. IEEE Trans. Evol. Comput. 2002, 6, 58−73. (27) Eberhart, R. C.; Kennedy, J. A New Optimizer Using Particle Swarm Theory. In Proceedings of the Sixth International Symposium on Micro Machine and Human Science; IEEE: New York, 1995; pp 39−43. (28) Sobol, I. M. The Distribution of Points in a Cube and the Approximate Evaluation of Integrals. USSR Comput. Math. Math. Phys. 1967, 7, 86−112. (29) Gentle, J. E. Random Number Generation and Monte Carlo Methods, 2nd ed.; Springer: New York, 2003. (30) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in FORTRAN 77: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992. (31) Matsumoto, L.; Nishimura, T. Mersenne Twister: A 623Dimensionally Equidistributed Uniform Pseudo-Random Number Generator. ACM Trans. Model. Comput. Simul. 1998, 8, 3−20. (32) Henderson, N.; Barufatti, N. E.; Sacco, W. F. The Least Dot Products Method: A New Numerical Paradigm for Phase Stability Analysis of Thermodynamic Mixtures. Chem. Eng. Sci. 2011, 66, 5684− 5702. (33) Shvarts, A. V.; Efremova, G. D. Higher-Order Critical Phenomena in the System Ethanol-Water-Carbon Dioxide. Russ. J. Phys. Chem. 1970, 44, 614−615.

4939


Prediction of Tricritical Points in Ternary Mixtures by Global

Recommend Documents