Calculation of Critical Points of Thermodynamic Mixtures with

1872

Ind. Eng. Chem. Res. 2010, 49, 1872–1882

Calculation of Critical Points of Thermodynamic Mixtures with Differential Evolution Algorithms Ne´lio Henderson,1,† Wagner F. Sacco,‡ Nelza E. Barufatti,† and M. M. Ali§ Thermodynamics and Optimization Group, Instituto Politećnico, UniVersidade do Estado do Rio de Janeiro, 28601-970, NoVa Friburgo-RJ, Brazil, UniVersidade Federal do Oeste Para´, 68040-070, Santare´m-PA, Brazil, and School of Computational and Applied Mathematics, Witwatersrand UniVersity, PriVate-Bag-3, Johannesburg, South Africa

In this article, we use differential evolution (DE) algorithms for the first time in the prediction of critical points of thermodynamic mixtures. Despite the name, DE algorithms are direct search populational techniques that do not use derivatives. We work with four different versions of DE algorithms and propose a variant for calculating more than one critical point. The differential evolution algorithms are tested on several multicomponent systems, and the new variant is capable of obtaining multiple solutions on different petroleum fluids. 1. Introduction A critical state of a thermodynamic mixture is an equilibrium state that is essentially between regions of stability and instability with regard to the number of phases. A critical state is characterized by its critical point, that is, by coordinates that describe the position of this state in the thermodynamic configuration space of the system. In liquid-vapor equilibria, for example, a critical point can be defined as the state at which vapor and liquid phases in equilibrium become indistinguishable. Unlike for pure substances, in mixtures, a critical point is not (necessarily) characterized by the highest pressure and temperature at which the system can exist in two phases. Thus, along a critical region, a peculiar phase transition occurs, called a second-order phase transition,1 which differs from the first-order transition that happens along the phase boundary, through bubble points and dew points. A second-order phase transition can engender complex thermodynamic phenomena, such as retrograde vaporization phenomenon, where the dew-point curve can exhibit complicated geometric forms when double-retrograde phenomena are present.2-5 To calculate critical points of multicomponent systems, several methods have been developed.6-16 In particular, the method of Hicks and Young8 is very interesting. This method is based on the idea that critical points are the intersections of some paths on the volume-temperature plane, the loci of the zeros of two criticality conditions. To locate these points, only one of the loci need be tracked. Thus, from initial points, the solutions are located by following the first condition locus and monitoring the sign of the second criticality condition. If the second condition changes sign, a solution has been detected. Hence, this solution can be located more precisely by an appropriate bisection technique. Although not mathematically guaranteed, as shown by Hicks and Young,8 in practice, this method can determine (with success) all critical points of binary mixtures. 1 To whom correspondence should be addressed. E-mail: [email protected]. † Universidade do Estado do Rio de Janeiro. ‡ Universidade Federal do Oeste Para´. § Witwatersrand University.

During the late 1970s and throughout the 1980s, critical-point thermodynamics saw great progresses. Peng and Robinson9 developed numerical strategies for predicting critical points from Gibbs criticality criteria based on the calculations of two determinants associated with the Gibbs energy. Later, Baker and Luks10 accomplished similar tasks using determinants associated with the Helmholtz energy. The evaluation of these determinants is a computationally expensive part of these two strategies. In fact, Baker and Luks10 observed that the computer time required for the solution of a critical-point problem increases exponentially with the number of components in the mixture. To decrease computation costs, Heidemann and Khalil11 obtained other forms for the criticality criteria using a Taylor expansion of the Helmholtz energy. They emphasized that this Taylor expansion of the Helmholtz energy is based on the thermodynamic concepts set forth by Gibbs.17 Following this path, criticality criteria based on the Taylor expansion of the tangent-plane distance function of the Gibbs energy was developed by Michelsen.13 This function (proposed originally by Gibbs and deduced formally by Baker et al.18) is the key to the modern stability analysis of phase equilibrium (see Michelsen19). From the Heidemann-Khalil formulation, a deterministic method that is guaranteed to find all critical points of a mixture was presented by Stradi et al.,16 who used an interval Newton/ generalized bisection algorithm. The methods of Hicks and Young8 and the interval Newton/ generalized bisection algorithm presented by Stradi et al.16 are frequently used. In fact, reliability is an inherent property of these algorithms, which motivates the use of both methodologies. However, as is well-known, these two methods are not computationally efficient for mixtures with large numbers of components. Henderson et al.20 developed a new methodology for the calculation of critical points of mixtures. They deduced critical conditions using a Taylor expansion from a slight modification of the Gibbs tangent-plane criterion and formulated the criticalpoint calculation, for the first time, as an optimization problem. In a subsequent work,21 they particularized the methodology for binary mixtures, illustrating its use in connection with the classification proposed by van Konynenburg and Scott.22 Such an approach offers many advantages, among them the formulation of critical-point problems from a two-variable objective

10.1021/ie900948z  2010 American Chemical Society Published on Web 01/11/2010

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010

function and the utilization of robust optimization algorithms that are easy to implement and do not require differentiability properties. Henderson et al.20,21 used the simulated annealing (SA) algorithm. This optimization method is very robust. However, for mixtures with large numbers of chemical species, the SA algorithm can be expensive. Several heuristic techniques for global optimization mimicking biological evolution have emerged in the literature. Among the most representative algorithms, we can mention the genetic algorithms23 (GAs), particle swarm optimization24 (PSO), and a new class of evolutionary methods called differential evolution (DE) algorithms.25,26 For different theoretical and practical problems, comparative studies have shown that the performance of DE algorithms is clearly better than that of the other two algorithms (GA and PSO); see, for example, Ali and To¨rn,27 Xu and Li,28 Pandura and Brizuela,29 and Sacco et al.30 To work with many chemical components, in this article, we use DE algorithms in the prediction of critical points. Despite the name, these robust and efficient algorithms are optimization techniques that do not use derivatives. A considerable number of variations of the differential evolution approach have been proposed in the literature.25-27,31-33 In chemical engineering, Wang and Jing34,35 introduced a hybrid method of differential evolution to obtain the maximum productivity in fermentation processes. Recently, Babu and Angira36 developed a modification of the original DE technique focused on the optimization of nonlinear chemical processes. Srinivas and Rangaiah37 used DE to solve phase equilibrium problems. Here, we use four different versions and propose a variant for calculating more than one critical point. DE algorithms are tested on several multicomponent systems, and the new variant is capable of obtaining multiple solutions on different petroleum fluids. All mixtures were modeled using the Peng-Robinson equation of state and the classical one-fluid van der Waals mixing rule.38 2. Optimization Problem 2.1. Formulation. Let I ) {x ) (x1, ..., xr-1) ∈ Rr-1; 0 < r-1 xi < 1, i ) 1, ..., r - 1 and ∑i)1 xi < 1} be a simplex in Rr-1, and let d: [0, +∞) × [0, +∞) × I f R, the modified tangentplane distance function, be given by

r-1

∇x2d(T, P, z)·u2 )

0 i

i

r

{

(1)

Here, µi ) µi(T,P,x) represents the chemical potential of component i () 1, ..., r) in mixture, and for all i, by definition, we have µi0 ) µi(T,P,z), where z is the vector of the (known) mole fractions of the chemical components in mixture. T is the temperature, and P the pressure. In agreement with Henderson et al.,20 given a mixture with composition z ∈ I, consider the functions q(T, P) ) ∇x d(T, P, z)·u (T, P, z)

(2)

c(T, P) ) ∇x3d(T, P, z)·u3(T, P, z)

(3)

2

where

2

(5)

(6)

To solve this nonlinear system, we should emphasize that, in each stage of any iterative method, the functions q(T,P) and c(T,P) are computed as indicated in algorithm 1. Algorithm 1 Step 1. Compute the Hessian matrix ∇x2d(T,P,z). Step 2. Compute u, an unitary eigenvector of ∇x2d(T,P,z) associated with the smallest eigenvalue of this Hessian matrix. Step 3. Set q(T,P) ) ∇x2d(T,P,z) · u2(T,P,z). Step 4. Compute ∇x3d(T,P,z) · u3(T,P,z). Step 5. Set c(T,P) ) ∇x3d(T,P,z) · u3(T,P,z) Note that, in step 2, the eigenvector calculation occurs on (r - 1) × (r - 1) matrices. Thus, for a binary mixture (r ) 2), the calculation is simplified. In this case, we have u ) 1. Hence, we obtain q(T,P) ) ∂2d(T,P,z)/∂x12 and c(T,P) ) ∂3d(T,P,z)/ ∂x13. In the present work, eigenvector calculations were performed using the classical Jacobi method,39 and the chemical potential derivatives were analytically obtained from Peng-Robinson equation of state with classical one-fluid van der Waals mixing rule. In agreement with Henderson et al.,20 the cubic form was approximated by central differences as

0 r

[µr(T, P, x) - µr0]

∂3d(T, P, z) uuu ∂xi∂xj∂xk i j k i,j,k)1

∑

q(T, P) ) 0 c(T, P) ) 0 s.t. (T, P) ∈ [Tmin, Tmax] × [Pmin, Pmax] ⊂ R2

r-1

i

(4)

j

are the second- and third-order tensors, respectively, of the derivatives of d(T,P,x) with regard to x, evaluated at the point (T,P,z) and applied in u, the unitary eigenvector of the Hessian matrix ∇x2d(T,P,z) associated with the smallest eigenvalue. Let [Tmin,Tmax] × [Pmin,Pmax] ⊂ R2 be a box where we can find a critical point (T,P) of a mixture with composition z. The critical-point calculation problem can be formulated as the following nonlinear system, with two equations and two variables

∑ x {[µ (T, P, x) - µ ] - [µ (T, P, x) - µ ]} + i)1

i j

i

r-1

∇x3d(T, P, z)·u3 )

∇x3d(T, P, z)·u3(T, P, z) ≈ d(T, P, x) )

2

P, z) uu ∑ ∂ d(T, ∂x ∂x

i,j)1

1873

1 [u·∇xd(T, P, z + δu) 6δ2 u·∇xd(T, P, z - δu)]

(7)

where δ > 0 is a sufficiently small scalar. We adopt δ ) 10-4|z|, where |z| is the Euclidian norm of z ) (z1, z2, ..., zr). The unitary vector u used in eq 7 is calculated in step 2 of algorithm 1. If the nonlinear system in eq 6 has a solution, then solving this system is equivalent to finding a global minimum of the objective function h(T,P) ) q(T,P)2 + c(T,P)2. Therefore, we implement the critical-point calculation as the following global optimization problem

{

given z ∈ I min h(T, P) ) q(T, P)2 + c(T, P)2 s.t. (T, P) ∈ [Tmin, Tmax] × [Pmin, Pmax] ⊂ R2

(8)

2.2. Why Global Optimization Techniques? In recent years, the field of global optimization has been developing rapidly. This progress has been an aid in solving many real-life

1874


Figure 1. Surface of the objective function for the mixture C2H6/n-C7H16/CH4.

Figure 2. Surface of the objective function for the mixture C2H6/n-C6H14/CH4.

problems that are difficult to optimize by traditional mathematical methods.40 Multicomponent mixtures can exhibit multiple critical points. Hence, the objective function h(T,P) described in eq 8 can present more than one global minimum. Moreover, in agreement with Henderson et al.,20 if the critical-point problem has only one solution, the optimization problem in eq 8 can also present several local minima, hindering the computation of the physically correct solution. In fact, as an illustration, Figure 1 presents the surface of the objective function h(T,P) for a ternary mixture containing C2H6 (component 1), n-C7H16 (component 2), and CH4 (component 3), with composition z1 ) 0.16392, z2 ) 0.0081372, and z3 ) 0.8279428. One can see that this function presents one global minimum and several local minima, all essentially in the same direction.

A more complex geometry is shown in Figure 2, where we consider a ternary mixture containing C2H6 (component 1), n-C6H14 (component 2), and CH4 (component 3), with z1 ) 0.10206, z2 ) 0.010982, and z3 ) 0.886958. In this example, the objective function h(T,P) has (seemingly) seven local minima, but only one global minimum. For effective visualization, in Figures 1 and 2, we consider h(T,P) ) 100, whenever h(T,P) > 100. However, we emphasize that the objective function h(T,P) assumes relatively high values in the vicinity of the critical point. There is no doubt that these optimization problems can be difficult to solve for many local optimization methods. Thus, such geometric difficulties justify the use of global optimization methods in the critical-point calculation of thermodynamic mixtures.


Moreover, many global optimization methods, such as differential evolution algorithms, are direct search methods. These optimization techniques avoid the calculation of derivatives of the objective function, thus diminishing the algebraic complexity of the optimization routine and becoming flexible for use with several thermodynamic models.

1875

(essentially) DE should change the population at each iteration, independently of the choice of the value of F. In this work, after intensive computational tests, we consider CR ) 0.4 and F ) 0.95. In addition, in step b of algorithm 2, we adopt ε ) 10-12.

3. Differential Evolution Algorithms 3.1. Original DE Algorithm. The differential evolution (DE) algorithm25,26 is an evolutionary direct search populational method, in which the optimality search is performed by mutating, crossing over, and selecting individuals along generations. Originally, the DE algorithm was designed for global optimization problems described as

{

min f ) f(x) s.t. x ∈ [a, b] ⊂ Rn

(9)

where [a,b] ) [a1,b1] × [a2,b2] × · · · × [an,bn] is a box in n-dimensional Euclidian space. All of the steps in DE are listed in algorithm 2 later in this section. First, an initial population (with np individuals) is randomly chosen. Then, inside a loop, the evolutionary process continues until a stopping criterion is satisfied. The stopping criterion considered in algorithm 2 (fbest e ε, where fbest is the best value of f) is particularly appropriate for critical-point calculations (see eq 6). Thus, once initialized, DE mutates and recombines the population to produce a population of np trial vectors. In the kth mutation (see step c2) a mutant vector is generated for each individual i as follows yi ) xr(k)1 + F· xr(k)2 - xr(k)3

[

]

(10)

where r1, r2, and r3 are random indexes that are mutually different from each other and different from index i and F is the scale factor, generally a contraction parameter in the range (0,1), that controls the rate at which the population evolves. The first random vector xr(k) is known as the base vector. 1 Therefore, in mutation, the base vector is altered by the addition of the weighted difference between the other vectors with indexes r2 and r3, which supplies the search direction. After mutation, the population goes through crossover (see step c3). In this stage, component j of trial vector zi is found from its parents xi(k) and yi according to the rule zi,j )

{

yi,j, if {rand(j)[0, 1] e CR} or j ) rand(i){1, ..., np} x(k) i,j , otherwise (11)

where rand(j)[0,1] is the jth evaluation of the random number generator in [0,1]; CR ∈ [0,1] is the crossover constant, which controls the fraction of parameter values that are copied from the trial vector; and rand(i){1, ..., np} is a randomly chosen index in {1, ..., np}. Note that the alternative j ) rand(i){1, ..., np} assures that at least one component will receive a mutated value. If zi ∉ [a,b], then we consider f(zi) ) fbig, where fbig > 0 is a sufficiently large constant. Finally, we have the selection process (step d in algorithm 2) that defines the population of the next generation. If vector zi yields a smaller objective function value than xi(k), then xi(k+1) is set to zi; otherwise, the old value xi(k) is retained. Note that the parameter F can be seen as the step length used to produce the mutants described by eq 10. On the other hand, the parameter CR determines the probability with which a mutant is accepted. Thus, for example, if CR ≈ 1, then

The critical-point problem is a two-dimensional optimization problem. Thus, the initial generation of the population can be built in a more appropriate particular way. Here, we divided the domain [Tmin, Tmax] × [Pmin, Pmax] ⊂ R2 into four quadrants of same size, and we generated np/4 individuals separately in each quadrant. Note that the population size (np) in this case should be chosen as multiple of four; we used np ) 20, 40, and 80. 3.2. Other Versions of the DE Algorithm. As originally proposed by Storn and Price,25 DE is a very robust algorithm, but for some problems, its convergence can be slow.32 To overcome this deficiency, several versions have been introduced in the literature. In this section, we describe three DE versions that are used here and compared to the original version. From algorithm 2, we note that, if f(zi) < f(x(k) i ), the offspring is not accepted immediately after its creation. Recently, based on this observation, Gong et al.33 suggested a modification whereby, if trial vector is better that xi(k), then zi is immediately copied into the current population. Thus, zi can be selected among the three solutions with indexes r1, r2, and r3 and contribute to the creation of better solutions. The variations proposed by Gong et al.33 for steps 3c and 3d of algorithm 2 are described in Figure 3. Here, algorithm 2 with these variations will be called the DE2 algorithm. To avoid the slowing of the convergence of the canonical DE algorithm as the region of the global minimum is approached, Kaelo and Ali32 incorporated a new idea regarding the mutation mechanism. They proposed a tournament selection to obtain the base vector, which can be mathematically described in the kth mutation as follows: Select uniform randomly distributed indexes r1, r2, and r3 in the set {1,2, ..., np} such that r1 * r2 * r3 * i. Let the subindex a ∈ {1, 2, 3} be such (k) (k) (k) that x(k) ra ) arg min{f(xr1 ), f(xr2 ), f(xr3 )}. Generate the offspring

1876


yi ) xr(k)a + F· xr(k)b - xr(k)c

[

]

(12)

where b, c ∈ {1, 2, 3} - {a} and b * c. (k) (k) Thus, given x(k) r1 , xr2 , and xr3 , the base vector is the best vector (k) among these three. In this new strategy, we note that p(k) i )(xrb - xr(k) ), the vector that defines the search direction of DE c algorithm, does not incorporate any information on x(k) ra , the best vector. To allow this eventuality, we generated the offspring simply as yi ) xr(k)a + F· xr(k)2 - xr(k)3

[

]

(13)

In agreement with Kaelo and Ali,32 algorithm 2 with this new mutation rule is referred to as the differential evolution algorithm with random localization (DERL). We also consider DERL2 algorithm, obtained from the DE2 algorithm together with this new mutation mechanism. The main differences among the four versions of the DE algorithm considered here are summarized in Table 1. 3.3. DE for Calculation of More Than One Critical Point. In this section, we introduce a DE variant designed for multimodal problems in the form shown in eq 9, considering the multimodal case, where a first solution x* ) (x*1 , ..., x*n ) ∈ [a, b] has been determined. To obtain another global minimum x**, we solved again the optimization problem with an alteration in the mutation step. Essentially, in this variant, the contraction factor F ∈ (0, 1) is replaced by an expansion factor E ∈ (1, +∞). At the ith (i ) 1,2, ..., np) mutation of the kth iteration of (k) DERL2, let x(k) s be the base vector and pi be the direction vector (k) (k) (k) in eq 13, that is, xs(k) ) xr(k) and p ) (x i r2 - xr3 ). The scaling a parameter of the mutation in eq 13 is determined by the value ξ, where (k) (k) (k) ξ ) min{|xs,1 - x*|, 1 |xs,2 - x*|, 2 ..., |xs,n - x*|} n

(14)

If ξ > σ, then eq 13 is given by yi ) xs(k) + F·p(k) i

(15)

otherwise

Table 1. Main Differences among the Four Versions of DE procedure

DE

DE2

DERL

DERL2

standard mutation mutation by tournament of Kaelo and Ali32 standard offspring acceptance modified offspring acceptance of Gong et al.33

yes no

yes no

no yes

no yes

yes

no

yes

no

no

yes

no

yes

on the size of the feasible region. For our case, we used E ) 100. The value of σ will be discussed in the next section. In the present article, DERL2-Mult denotes the DERL2 algorithm with this variation in the mutation rule. 4. Results and Discussion All numerical results in this work were obtained using an AMD Athlon 64 X2 Dual Core Processor 400+ PC with 2 Gb of RAM. The DE algorithm and its versions were implemented in the FORTRAN 77 language, with the Ifort-Intel compiler for Linux. For the stochastic part of the differential evolution algorithms, we used the pseudorandom-number-generating algorithm developed by Matsumoto and Nishimura.41 4.1. Comparing the Different Versions. In Table 2, we consider 19 mixtures, n-alkane and N2-n-alkane systems, for which experimental critical points were available in the literature. Note that these mixtures include 11 three-component, four four-component, three five-component, and one six-component mixture. For all mixtures in Table 2, the input data for component i () 1, ..., r), namely,the critical temperature (Tci), critical pressure (Pci), and acentric factors (ωi), were taken from Reid et al.42 In the mixing rule, all binary interaction parameters were considered to be null. For these mixtures, the optimization problems were formulated on the domain [Tmin, Tmax] × [Pmin, Pmax], whose upper and lower limits were defined as follows r

Tmin )

∑ i)1

yi ) xs(k) + E·p(k) i

(16)

Here, σ > 0 is the neighborhood radius that defines the expansion process. If two solutions x*, x** ∈ [a, b] have already been obtained, then in the determination of the third minimum x*** ∈ [a, b], the distance ξ is calculated as ξ ) min{ξ1,ξ2}, where (k) (k) (k) ξ1 ) min{|xs,1 - x*1 |, ..., |xs,n - x*n |} and ξ2 ) min{|xs,1 - x** 1 |, (k) ** ..., |xs,n - xn |}, and so on. If E is relatively large, then it is possible that zi ∉ [a,b]. In this case, we will have f(zi) ) fbig, and zi will certainly not be selected. The value of E is dependent

Figure 3. Modifications in the DE algorithm by Gong et al.33

r

ziTci, Tmax ) 4.0

∑zT

r

Pmin ) 0.4

∑ i)1

i ci

i)1

(17)

r

ziPci,

Pmax ) 4.0

∑zP i)1

i ci

(18)

To solve these critical-point problems, we used the original DE algorithm as well as the DE2, DERL, and DERL2 versions. Table 3 lists the experimental critical points and the calculated values obtained using the evolution differential algorithms, from the Peng-Robinson equation of state.38 We emphasize that the numerical simulations in this section do not have the objective of validating the thermodynamic model. In fact, this was already done by Peng and Robinson.9 Our main objective was to compare the different versions of the differential evolution algorithm in the context of criticalpoint calculations, with well-known mixtures taken from the literature. Thus, the experimental results (see Table 3) are included here, essentially, to allow a judgment of the present numerical implementation, including the routine of the classical Peng-Robinson model. As one can see, the numerical results are in agreement with the experimental values. The performances of the DE algorithms are compared in Table 4, where the number of objective function evaluations (NFE) and the elapsed time (in seconds) are included. We emphasize that all of the differential evolution algorithms


1877

Table 2. Compositions of 19 Multicomponent Systems with n-Alkanes and Nitrogen mixture no.

no. of components

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

3 3 3 3 3 3 3 4 4 4 5 3 3 3 3 4 5 5 6

C1

C2

C3

0.4290 0.7260 0.5140 0.8010 0.6120 0.6150 0.3414 0.6168 0.2542 0.2019 0.415 0.360 0.4530 0.4115 0.4345 0.6626 0.7057 0.1015

0.2029

0.0835 0.1093 0.0669 0.3573

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

calculated Tc(K)

0.3421

440.24 6295.14 389.74 7485.63 404.41 6232.18 395.87 8306.74 425.94 7050.42 420.67 6946.50 404.24 5541.39 425.47 7374.14 410.60 5056.71 419.66 4422.20 390.37 6860.06 327.68 8527.26 328.17 8985.02 321.58 8983.00 321.93 9347.98 316.27 8874.06 321.64 14489.9 320.31 14574.1 381.26 6440.71

Tc (K) Pc (kPa) 438.15 385.92 400.37 391.48 421.48 415.92 397.15 423.15 405.87 417.92 387.03 322.03 322.70 313.70 313.70 313.70 310.53 308.42 376.42

6612 7605 6405 8101 7156 7060 5602 7412 5113 4506 7220 8674 9204 9232 9797 8963 13748 13700 6536

0.1376 0.2554 0.3316 0.2038

0.2547 0.4858 0.2033 0.542 0.545 0.5005 0.5030 0.4330 0.1057 0.0413 0.2629

experimental

Pc (kPa)

n-C5

n-C6

0.373 0.171 0.412

Table 3. Calculated Critical Points and Corresponding Experimental Data mixture no.

n-C4

ref Ekiner and Thodos43 Ekiner and Thodos43 Ekiner and Thodos43 Ekiner and Thodos44 Ekiner and Thodos44 Ekiner and Thodos44 Etter and Kay45 Ekiner46 Etter and Kay45 Etter and Kay45 Etter and Kay45 Yarborough and Smith47 Yarborough and Smith47 Yarborough and Smith47 Yarborough and Smith47 Yarborough and Smith47 Hanson and Brown48 Hanson and Brown48 Etter and Kay45

obtained the same numerical solutions reported in Table 3, using the same values of ε ) 10-12 (see step 2 in algorithm 2) and np ) 40.

0.0640 0.2710 0.2960 0.3165 0.0726 0.2357 0.1213 0.1881

n-C7

N2

0.1980 0.1030 0.0740 0.1350 0.1170 0.0890 0.1730 0.0613 0.0430 0.0950 0.0465 0.0855 0.0490

0.0508 0.1794

0.0616 0.1353 0.0657

0.0608 0.0332

Plots with the relative errors taken with respect to the experimental data are shown in Figure 4. In most of the calculated results, we note that the relative error is smaller that 2% for the critical temperature and smaller than 3% for the critical pressure. For the systems studied, Figure 5 clarifies that, essentially, DE, DE2, DERL, and DERL2 do not exceed NFE ) 6000, 5000, 4500, and NFE ) 2500, respectively. The DERL2 algorithm is the most efficient one. In fact, this version did not exceed 0.02 s for all considered mixtures, and in various cases, we obtained some processing times less than or equal to 0.01 s. In the present work, all calculated critical points were tested for global phase stability, using the robust methodology presented by Henderson et al.49 4.2. Comparison with the Results of Stradi et al. Mixtures 1-6 and 12-15 in Table 2 were studied by Stradi et al.16 Consequently, in this section, we compare the critical points of these mixtures, as calculated by Stradi et al.,16 with the values obtained by the DE algorithms. For this purpose, we employed the nonzero binary interaction parameters (ki,j) shown in Table 5, the same parameters as used by Stradi et al.16

Table 4. Performances of the Differential Evolution Algorithms DE

DE2

DERL

DERL2

mixture no.

no. of components

NFE

time (s)

NFE

time (s)

NFE

time (s)

NFE

time (s)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

3 3 3 3 3 3 3 4 4 4 5 3 3 3 3 4 5 5 6

5880 5640 5560 5240 5480 5880 12680 5640 6040 2360 5960 5480 5080 4520 5080 5480 5160 4680 6040

0.02 0.02 0.02 0.02 0.02 0.02 0.04 0.02 0.02 0.01 0.04 0.02 0.01 0.01 0.02 0.02 0.03 0.03 0.05

5160 4200 4600 4200 3720 4440 10760 4600 6680 5480 4600 4120 3720 4200 3320 3880 4360 3320 5720

0.02 0.01 0.01 0.01 0.01 0.01 0.04 0.02 0.03 0.02 0.03 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.05

5320 4120 3480 3720 3400 3400 7640 3400 4360 2360 3720 3240 5680 4520 2840 3160 2760 2600 3720

0.02 0.01 0.01 0.01 0.01 0.01 0.03 0.01 0.02 0.01 0.03 0.01 0.02 0.01 0.01 0.01 0.02 0.02 0.03

3960 2120 2040 3320 2120 1960 6200 2280 2520 2680 2440 1720 2200 1960 1720 1880 1560 1560 2440

0.01 0.00 0.00 0.01 0.00 0.00 0.02 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.02

1878


Figure 4. Relative errors for critical temperature and critical pressure.

Figure 5. Numbers of objective function evaluations for the different versions of DE. Table 5. Binary Interation Parameters (ki,j) from Stradi et al.16

N2 C1 C2 C3 n-C4 n-C5 n-C7

N2

C1

C2

C3

n-C4

n-C5

n-C7

0.0000 0.0311 0.0000 0.0852 0.0000 0.0000 0.0000

0.000 0.000 0.0140 0.000 0.000 0.000

0.0000 0.0000 0.0096 0.0078 0.0067

0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0033

0.0000 0.0074

0.0000

Table 6 reports the critical points obtained in this work and those obtained by Stradi et al.16 To facilitate the comparison, we converted the pressure values to the same units as used by these authors (bar instead of kPa). Observing Table 6, we can

affirm that the two methodologies obtained essentially the same critical points for all of the systems considered. Beyond this comparison, we can note the sensitivity of the critical points with respect to ki,j. In fact, the results in Tables

Ind. Eng. Chem. Res., Vol. 49, No. 4, 2010 Table 6. Comparison with the Results of Stradi et al. this work

16

Table 7. Binary Mixture (CH4/H2S)

Stradi et al.16

mixture no.

Tc (K)

Pc (bar)

Tc (K)

Pc (bar)

1 2 3 4 5 6 12 13 14 15

439.54 388.64 403.53 394.71 424.76 419.54 327.46 328.15 321.26 321.75

63.2 74.8 62.5 82.9 70.4 69.38 87.0 92.5 91.8 96.2

439.54 388.64 403.54 394.71 424.77 419.54 327.47 328.15 321.27 312.76

63.2 74.9 62.5 82.9 70.4 69.39 87.0 92.5 91.8 96.2

3 and 6 clearly indicate that the calculated points are sensitive to the chosen values for the binary interaction parameters. 4.3. Computing More Than One Critical Point. In the considered ranges, none of the mixtures studied so far in this work present more than one critical point. On the other hand, in the next subsections, all mixtures have two or more critical points. In these cases, we will use the variation of the differential evolution algorithm proposed in section 3.3. The robustness and efficiency of the DERL2 algorithm indicate the choice of the DERL2-Mult algorithm (see section 3.3) to solve such multimodal problems, a difficult task for any optimization method.20 In the next problems, with the exception of the samples in subsection 4.3.1, the upper and lower limits are defined as Tmin ) 100 K, Pmin ) 1000 kPa,

1879

Tmax ) 650 K

(19)

Pmax ) 60000 kPa

(20)

4.3.1. Binary Mixture. Here, calculations were performed on the binary mixture CH4/H2S with different compositions. Freitas et al.21 calculated the critical points of that mixture using the simulated annealing algorithm. The critical temperature, critical pressure, and acentric factors were taken from Reid et al.,42 and we used kCH4,H2S ) 0.08. In this subsection, we considered the following ranges: [Tmin, Tmax] ) [200, 500], [Pmin, Pmax] ) [400, 25000] if zCH4 ∈ [0.49, 0.52] and [Pmin, Pmax] ) [600, 150000] otherwise, where the temperatures are given in Kelvin and the pressures are in kilopascals. Table 7 reports the results computed by DERL2-Mult. As observed before,21 this mixture can present two and up to three critical points, which are generally stable. The DERL2-Mult algorithm did not exceed the value NFE ) 5080 in approximately 0.003 s, and for many samples, we obtained NFE ≈ 3000. Using the simulated annealing algorithm, Freitas et al.21 obtained the same results, with NFE ≈ 9500. 4.3.2. Eight-Component Mixture of Nagarajan et al. We also considered the eight-component mixture studied initially by Nagarajan et al.50 The composition and the equation-of-state properties for this reservoir fluid are given in Table 8, including the nonzero binary interaction parameters (ki,j). Table 8 was obtained from Saber and Shaw.51 They did not calculate critical points, but they did mention the existence of one critical point near T ) 353 K and P ) 38500 kPa. In fact, this point was used by Saber and Shaw51 to perform phase equilibrium calculations near a critical point. The results obtained with the DELR2-Mult algorithm using σ ) 8 are given in Table 9. As shown in this table, we found two stable critical points, with the one at Tc ) 362 K and P ) 39878.61 kPa probably being the point mentioned by Saber and Shaw.51 We note that the DERL2-Mult algorithm did not exceed the value NFE ) 3560 in 0.062 s.

CH4 (%)

multiple CPs

Tc (K)

Pc (kPa)

52.0

first CP (unstable) second CP (stable) first CP (stable) second CP (stable) first CP (stable) second CP (stable) first CP (stable) second CP (unstable) third CP (stable) first CP (stable) second CP (stable) third CP (stable) first CP (stable) second CP (stable) third CP (stable) first CP (stable) second CP (stable) third CP (stable) first CP (stable) second CP (stable) third CP (stable) first CP (stable) second CP (stable) third CP (stable)

254.75 269.54 244.72 277.58 227.65 287.06 224.38 205.35 288.22 223.15 205.88 288.60 221.81 206.53 288.98 220.30 207.37 289.35 208.51 218.50 289.72 210.38 215.99 290.80

15 096.49 14 664.76 16 275.68 14 606.88 24 113.77 14 520.76 27 583.49 148 027.17 14 505.38 29 199.99 132 020.69 14 500.33 31 206.92 116 929.58 14 494.74 33 822.62 102 381.74 14 489.36 87 722.65 37 546.55 14 483.75 70 971.89 44 236.01 14 478.10

51.0 49.0 48.7 48.6 48.5 48.4 48.3 48.2

NFE np 3240 5080 2600 3880 2520 5080 2840 2920 5080 3000 3480 5080 2520 3080 5080 2600 2760 5080 3240 3400 5080 2600 3160 5080

σ 5.0 11.0 30.0 1.0 46.0 1.0 45.0 2.0 4.5 1.0 45.0 1.0 45.0 1.0 45.0

40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40

Table 8. Eight-Component Mixture from Saber and Shaw51 component mole fraction Tci (K) Pci (kPa) C1 C2 C3 n-C4 n-C5 n-C6 C7-C16 C17+

0.6883 0.0914 0.0460 0.0333 0.0139 0.0152 0.0896 0.0223

190.60 305.40 369.80 425.20 469.60 507.40 606.28 825.67

4599 4883 4244 3799 3373 2968 2575.7 1458

ki,C7-C16 ki,C17+

ωi 0.008 0.0908 0.1520 0.1930 0.2510 0.2960 0.4019 0.7987

0.05 0.04 0.01 0.00 0.00 0.00 0.00 0.00

0.090 0.055 0.010 0.000 0.000 0.000 0.000 0.00

Table 9. DELR2-Mult Results for the Eight-Component Mixture multiple CPs

Tc (K)

Pc (kPa)

NFE

time (s)

np

σ

first CP (stable) second CP (stable)

259.91 362.65

42 477.52 39 878.61

1700 3560

0.03 0.62

20 40

8.0

Table 10. Thirteen-Component Mixtures of Ghorayeb et al.52 composition (%) component C30+ C25-C29 C20-C24 C16-C19 C13-C14-C15 C10-C11-C12 C7-C8-C9 C6 n-C5-i-C5 n-C4-i-C4 C3 C2 C1-CO2-N2

G1

G2

G3

0.49 1.34 0.59 0.74 0.88 1.26 2.97 0.65 0.93 2.28 3.78 8.15 76.96

0.70 1.91 2.03 1.82 1.72 2.02 4.04 0.80 1.08 2.48 3.89 8.06 69.45

0.40 3.89 3.34 2.60 2.23 2.42 4.49 0.86 1.11 2.47 3.82 7.84 64.53

Tci (K) Pci (kPa) 830.00 818.92 778.92 726.96 679.05 622.29 554.13 506.35 469.60 425.18 369.80 305.40 190.60

686 1108 1328 1679 2031 2453 2969 3396 3330 3799 4190 4883 4600

ωi

ki,C1-CO2-N2

1.505 1.399 1.168 0.940 0.776 0.611 0.428 0.299 0.251 0.193 0.152 0.098 0.008

0.024 0.020 0.020 0.019 0.015 0.010 0.009 0.000 0.000 0.000 0.000 0.000 0.000

Table 11. DELR2-Mult Algorithm Results for the 13-Component Mixture of Ghorayeb et al.52 mixture

multiple CPs

G1

first CP (unstable) second CP (stable) one CP (stable) one CP (stable)

G2 G3

Tc (K) Pc (kPa) NFE time (s) np 226.22 332.23 549.28 619.00

27 194.19 40 157.62 33 068.84 25 368.32

4522 4522 4522 4522

0.23 0.23 0.23 0.23

40 40 40 40

σ 3.0 -

4.3.3. Thirteen-Component Mixtures of Ghorayeb et al. Table 10 shows the composition, equation-of-state properties, and nonzero binary interaction parameters for the 13-component

1880


Table 12. Sixteen-Component Mixture with Carbon Dioxide and Pseudocomponents (PCs) component

mole fraction

CO2 N2 C1 C2 C3 i-C4 n-C4 i-C5 n-C5 C6 PC1 PC2 PC3 PC4 PC5 PC6

7.28017 1.32016 4.48307 1.10839 8.16854 9.90126 9.04865 4.34555 5.91325 9.13116 4.98710 4.43876 3.44657 2.62408 1.59272 6.25869

× × × × × × × × × × × × × × × ×

-1

10 10-3 10-2 10-2 10-3 10-4 10-3 10-3 10-3 10-3 10-2 10-2 10-2 10-2 10-2 10-3

Tci (K)

Pci (kPa)

ωi

ki,CO2

ki,N2

ki,C1

304.20 126.20 190.60 305.40 369.80 408.10 425.20 468.77 469.60 511.98 566.55 647.06 719.44 784.93 846.33 919.39

7380 3390 4600 4880 4250 3650 3800 3558 3370 3319 2530 1910 1420 1050 750 476

0.2250 0.0400 0.0115 0.0908 0.1454 0.1760 0.1928 0.2202 0.2273 0.2606 0.3884 0.5289 0.6911 0.8782 1.1009 1.4478

0.00 -0.020 0.075 0.080 0.080 0.085 0.085 0.085 0.085 0.095 0.095 0.095 0.095 0.095 0.095 0.095

0.00 0.00 0.08 0.07 0.07 0.06 0.06 0.06 0.06 0.05 0.10 0.12 0.12 0.12 0.12 0.12

0.000 0.000 0.000 0.003 0.010 0.018 0.018 0.025 0.026 0.036 0.049 0.073 0.098 0.124 0.149 0.181

Table 13. DELR2-Mult Algorithm Results for the 16-Component Mixture multiple CPs

Tc (K)

Pc (kPa)

NFE

time (s)

np

σ

first CP (stable) second CP (unstable)

590.69 292.49

20991.73 8899.79

5720 11920

0.61 1.24

40 80

40

mixtures of Ghorayeb et al.,52 here called mixtures G1, G2, and G3. According to Ghorayeb et al., these systems are unusual fluids of a gas-condensate field. The results calculated with the DELR2-Mult algorithm using np ) 40 and σ ) 3 are listed in Table 11. For mixture G1, we obtained two critical points, but one was unstable. Using an efficient numerical strategy, Hoteit et al.53 calculated only one critical point for all mixtures. Other than the above-mentioned unstable critical point, we observed that the results obtained with DELR2-Mult agreed with those of Hoteit et al.53 For these 13-component mixtures, each calculation with DELR2-Mult did not exceed 0.3 s. 4.3.4. Sixteen-Component Mixture with an Elevated Proportion of Carbon Dioxide. The composition in Table 12

was obtained from oil B of Shelton and Yarborough54 by CO2 addition. In this example, there is a high proportion of carbon dioxide (∼73%). The equation-of-state properties and the nonzero binary interaction parameters were obtained from Nichita et al.55 Two critical points calculated with the DELR2Mult algorithm are listed in Table 13, where one critical point is unstable. For this 16-component mixture, the first critical point was obtained in 0.61 s, using np ) 40. In this example, we note that the second calculation was a more difficult task, accomplished with np ) 80, NFE ) 11920, and time ) 1.24 s. 4.3.5. Reservoir Oil with Twenty-Nine Components. Table 14 lists data for the reservoir oil with 29 components studied by Nichita et al.56 In Table 14, ni represents the number of moles of component i in the mixture. Nichita et al.56 used this mixture to determine the phase envelope, the spinodal, the convergence locus, and the stability test limit locus in the presence of CO2 injection. The results calculated with the DELR2-Mult algorithm using σ ) 3.5 are reported in Table 15. We obtained four critical points. Note that three of the calculated critical points are

Table 14. Reservoir Oil Described by 29 Components component

ni

Tci (K)

Pci (kPa)

ωi

ki,N2

ki,CO2

ki,H2S

ki,COS

ki,C1

ki,C2

ki,C3

N2 CO2 H2S CH3S C2H6S COS C1 C2 C3 i-C4 n-C4 i-C5 n-C5 C6 C6H6 C7 C7H8 C8 C8H10(XYL) C8H10(EB) C9 C10 C11 C12 C13 C14 PC1 PC2 PC3

1.150 4.120 14.69 0.012 0.016 0.002 48.21 7.320 4.430 0.860 1.930 0.890 0.880 1.270 0.049 1.591 0.178 1.702 0.343 0.057 1.250 1.330 1.050 0.840 0.780 0.620 2.388 1.762 0.310

126.20 304.20 373.20 469.95 499.00 378.80 190.60 305.40 369.80 408.10 425.20 468.77 469.60 511.98 562.20 549.75 591.80 574.715 612.13 617.10 618.15 638.15 658.15 676.15 690.15 708.15 741.78 906.15 1113.15

3390 7380 8937 7230 5490 6350 4600 4880 4250 3650 3800 3558 3370 3319 4890 3063 4100 2830 3600 3610 2758 2585 2457 2320 2196 2078 1819 1680 1435

0.0400 0.2250 0.1000 0.1530 0.1910 0.1050 0.0115 0.0908 0.1454 0.1760 0.1928 0.2202 0.2273 0.2606 0.2120 0.2807 0.2630 0.3270 0.3230 0.3011 0.3440 0.3780 0.3960 0.4140 0.4320 0.4590 0.5129 0.6950 0.8300

0.00 -0.02 0.18 0.00 0.00 0.00 0.04 0.05 0.08 0.10 0.09 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10

-0.02 0.00 0.10 0.00 0.00 0.00 0.10 0.13 0.13 0.13 0.13 0.12 0.12 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10

0.180 0.100 0.000 0.000 0.000 0.000 0.000 0.150 0.090 0.075 0.050 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060 0.060

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.15 0.15

0.04 0.10 0.00 0.15 0.15 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.07 0.00 0.07 0.00 0.07 0.07 0.00 0.00 0.00 0.00 0.00 0.00 0.02 0.08 0.08

0.05 0.13 0.15 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.00 0.03 0.00 0.03 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.04

0.08 0.13 0.09 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.03 0.03


1881

Table 15. DELR2-Mult Algorithm Results for the 29-Component Mixture multiple CPs

Tc (K)

Pc (kPa)

NFE

time (s)

np

σ

first CP (unstable) second CP (unstable) third CP (unstable) fourth CP (stable)

230.60 235.05 240.18 424.87

3042.93 4214.52 5127.10 33941.51

2460 2100 3660 5560

0.90 0.74 1.21 1.86

20 20 20 40

3.5 3.5 3.5

unstable. The only stable critical point is located at T ) 424.87 K and P ) 33941.51 kPa. Nichita et al.56 reported that a critical point at T ) 412.89 K and P ) 33881.00 kPa. Thus, regarding the solutions indicated by these authors, the stable critical point calculated in the present work has relative errors equal to 0.18% for pressure and 2.9% for temperature. If this is caused by the approximation of using finite differences shown in eq 7, then this type of error can probably be decreased by using analytical expressions for the cubic-form calculation, mainly for mixtures with many components. All critical points in Table 15 were calculated in less than 2 s, corroborating the efficiency and robustness of the DELR2Mult algorithm. 5. Conclusions In this work, the calculation of critical points is described as a global optimization problem. As emphasized by Henderson et al.,20 this approach offers many advantages, including the formation of the critical-point problem as a two-variable objective function and the ease of implementation of robust direct search global optimization algorithms. Here, the optimization problems were solved using four different versions of the differential evolution algorithm: DE, DE2, DELR, and DELR2. In addition, we proposed a new variant for calculating more than one critical point, denoted DELR2-Mult. The differential evolution algorithms were tested with success on several petroleum fluids. In general, using the DELR2 algorithm, critical points could be calculated in a few hundredths of a second. Acknowledgment N.H. gratefully acknowledges the research grant provided by CNPq (Conselho Nacional de Desenvolvimento Cientı´fico e Tecnolo´gico, Ministry of Science & Technology, Brazil). The research by N.H. was carried out within the framework of Project PROCIENCIA-UERJ financed by FAPERJ. W.F.S. gratefully acknowledges the research grant provided by CNPq and FAPESPA (Fundaçaõ de Amparo a` Pesquisa do Estado do Para´, State of Para´, Brazil). The authors are grateful to the two anonymous reviewers for their suggestions and remarks that were decisive in improving this article. We are also grateful to Dr. Leá Freitas for the production of Figures 1 and 2. Literature Cited (1) Callen, H. B. Thermodynamics and an Introduction to Thermostatistics, 2nd ed.; Wiley: New York, 1985. (2) Raeissi, S.; Peters, C. J. On the Phenomenon of Double Retrograde Vaporization: Multi-Dew Point Behavior in the Binary System Ethane + Limonene. Fluid Phase Equilib. 2001, 191, 33. (3) Raeissi, S.; Peters, C. J. Simulation of Double Retrograde Vaporization Using the Peng-Robinson Equation of State. J. Chem. Thermodyn. 2003, 35, 573. (4) Raeissi, S.; Peters, C. J. Application of Double Retrograde Vaporization as an Optimizing Factor in Supercritical Fluid Separations. J. Supercrit. Fluids 2005, 115, 33. (5) Henderson, N.; Sacco, W. F. Prediction of Double Retrograde Vaporization by Hybrid Global-Local Optimization Factor Using Fuzzy Clustering Means. Chemical Product and Process Modeling 2008, 3, article 47.

(6) Spear, R. R.; Robinson, R. L., Jr.; Chao, K. C. Critical States of Ternary Mixtures and Equations of States. Ind. Eng. Chem. Fundam. 1971, 10, 588. (7) Hicks, C. P.; Young, C. L. The Gas-Liquid Critical Properties of Binary Mixtures. Chem. ReV. 1975, 75, 119. (8) Hicks, C. P.; Young, C. L. Theoretical Predictions of Phase Behavior at High Temperatures and Pressure for Non-Polar Mixtures: 1. Computer Solution Techniques and Stability Testes. J. Chem. Soc., Faraday Trans. 2 1977, 73, 597. (9) Peng, D.; Robinson, D. B. A Rigorous Method for Predicting the Critical Properties of Multicomponent Systems from an Equation of State. AIChE J. 1977, 23, 137. (10) Baker, L. E.; Luks, K. D. Critical Point and Saturation Pressure Calculations for Multipoint Systems. Soc. Pet. Eng, J. 1980, 20, 15. (11) Heidemann, R. A.; Khalil, A. M. The Calculation of Critical Points. AIChE J. 1980, 26, 769. (12) Michelsen, M. L. Calculation of Phase Envelopes and Critical Points for Multicomponent Mixtures. Fluid Phase Equilib. 1980, 4, 1. (13) Michelsen, M. L. Calculation of Critical Points and Phase Boundaries in the Critical Region. Fluid Phase Equilib. 1984, 16, 57. (14) Michelsen, M. L.; Heidemann, R. A. Calculation of Critical Points from Cubic Two-Constant Equations of State. AIChE J. 1981, 27, 521. (15) Wang, M. C.; Wong, D. S.; Chen, H.; Yah, W.; Guo, T. M. Homotopy Continuation Method for Calculation of Critical Loci of Binary Mixtures. Chem. Eng. Sci. 1999, 54, 3873. (16) Stradi, B. A.; Brennecke, J. F.; Kohn, J. P; Stadtherr, M. A. Reliable Computation of Mixture Critical Points. AIChE J. 2001, 47, 212. (17) Gibbs, J. W. On the Equilibrium of Heterogeneous Substances. Collected Works (October 1876-May 1877); Yale University Press: New Haven, CT, 1928; Vol. 1, p 55. (18) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng, J. 1982, 22, 731. (19) Michelsen, M. L. The Isotherm Flash Problem. Part 1. Stability Analysis. Fluid Phase Equilib. 1982, 9, 1. (20) Henderson, N.; Freitas, L.; Platt, G. M. Prediction of Critical Points: A New Methodology Using Global Optimization. AIChE J. 2004, 50, 1300. (21) Freitas, L.; Henderson, N.; Platt, G. M. Novel Approach for the Calculation of Critical Points in Binary Mixtures Using Global Optimization. Fluid Phase Equilib. 2004, 225, 29. (22) van Konynenburg, P. H.; Scott, R. L. Critical Lines and Phase Equilibria in Binary van der Waals Mixtures. Philos. Trans. R. Soc. London A: Phys. Sci. 1980, 298, 495. (23) Golberg, D. E. Genetic Algorithms in Search, Optimization, and Machine Learning; Addison- Wesley: Boston, 1989. (24) Kennedy, J.; Eberhart, R. Particle Swarm Optimization. Proc. IEEE Int. Conf. Neural Networks 1995, 4, 1942. (25) Storn, R. M.; Price, K. V. Differential EvolutionsA Simple and Efficient Heuristic for Global Optimization. J. Global Optim. 1997, 11, 341. (26) Price, K. V.; Storn, R. M.; Lampinen, A. J. Differential EVolution: A Practical Approach to Global Optimization; Springer: Berlin, 2005. (27) Ali, M. M.; Torn, A. Population Set Based Global Optimization Algorithms: Some Modifications and Numerical Studies. Comput. Oper. Res. 2004, 31, 1703. (28) Xu, X., Li, Y. Comparison between Particle Swarm Optimization, Differential Evolution and Multi-Parents Crossover. In Proceedings of the IEEE International Conference on Computational Intelligence and Security; IEEE Press: Piscataway, NJ, 2007; pp 124-127; doi 10.1109/cis.2007.37. (29) Panduco, M. A.; Brizuela, C. A. A Comparison of Genetic Algorithms, Particle Swarm Optimization and the Differential Evolution Method for the Design of Scannable Circular Antenna Arrays. Prog. Electromagn. Res. B 2009, 13, 171. (30) Sacco, W. F.; Henderson, N.; Rios-Coelho, A. C.; Ali, M. M.; Pereira, C. M. N. A. Differential Evolution Algorithms Applied to Nuclear Reactor Core Design. Ann. Nucl. Energy 2009, 36, 1093. (31) Fan, H. Y.; Lampinen, A. J. A Trigonometric Mutation Operation to Differential Evolution. J. Global Optim. 2003, 27, 105. (32) Kaelo, P.; Ali, M. M. A Numerical Study of Some Modified Differential Evolution Algorithms. Eur. J. Oper. Res. 2006, 169, 1176.

1882


(33) Gong, W.; Cai, Z.; Jiang, L. Enhancing the Performance of Differential Evolution Using Orthogonal Design Method. Appl. Math. Comput. 2008, 206, 56. (34) Wang, F. S.; Jing, C. H. Fuzzy-Decision-Making Problems of Fuel Ethanol Production Using a Genetically Engineered Yeast. Ind. Eng. Chem. Res. 1998, 37, 3434. (35) Wang, F. S.; Jing, C. H. Performance Analysis and Fuzzy Optimization of a Two-Stage Fermentation Process with Cell Recycling Including an Extractor for Lactic Acid Production. Chem. Eng. Sci. 2003, 58, 3753. (36) Babu, B. V.; Angira, R. Modified Differential Evolution (MDE) for Optimization of Non-linear Chemical Processes. Comput. Chem. Eng. 2006, 30, 989. (37) Srinivas, M.; Rangaiah, G. P. A Study of Differential Evolution and Tabu Search for Benchmark Phase Equilibrium and Phase Stability Problems. Comput. Chem. Eng. 2007, 31, 760. (38) Peng, D.; Robinson, B. D. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59. (39) Sewell, G. Computational Methods of Linear Algebra; Ellis Horwood: London, 1990. (40) Floudas, C. A.; Gounaris, C. E. A Review of Recent Advances in Global Optimization. J. Global Optim. 2009, 45, 3. (41) Matsumoto, L.; Nishimura, T. Mersenne Twister: A 623-dimensionally Equidistributed Uniform Pseudo-Random Number Generator. ACM Trans. Model. Comput. Simul. 1998, 8, 3. (42) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids; McGraw-Hill: New York, 1988. (43) Ekiner, O.; Thodos, G. Critical Temperatures and Pressures of the Ethane-n-Butane-n-Heptane System. J. Chem. Eng. Data 1968, 13, 304. (44) Ekiner, O.; Thodos, G. Critical Temperatures and Pressures of the Ethane-n-Pentane-n-Heptane System. J. Chem. Eng. Data 1966, 11, 457. (45) Etter, D. O.; Kay, G. Critical Properties of Mixtures of Normal Paraffin Hydrocarbons. J. Chem. Eng. Data 1961, 6, 409. (46) Ekiner, O. Critical Temperature and Pressure for Multicomponet Hydrocarbon Systems. Ph.D. Dissertation, Northwestern University, Evanston, IL, 1965.

(47) Yarborough, L.; Smith, L. R. Solvent and Driving Gas Compositions for Miscible Slug Displacement. Soc. Petr. Eng. J. 1970, 10, 298. (48) Hanson, G. H.; Brown, C. G. Vapor-Liquid Equilibrium in Mixtures of Volatile Paraffins. Ind. Eng. Chem. 1945, 37, 821. (49) Henderson, N.; de Oliveira, J. R., Jr.; Souto, H. A.; Marques, R. P. Modelling and Analysis of the Isothermal Flash and Its Calculation with the Simulated Annealing Algorithm. Ind. Eng. Chem. Res. 2001, 40, 6028. (50) Nagarajan, N. R.; Cullick, A. S.; Griewank, A. New Strategy for Phase Equilibrium and Critical Point Calculations by Thermodynamic Energy Analysis Part 1. Stability Analysis and Flash. Fluid Phase Equlib. 1991, 64, 191. (51) Saber, N.; Shaw, J. M. Rapid and Robust Phase Behaviour Stability Analysis Using Global Optimization. Fluid Phase Equlib. 2008, 264, 137. (52) Ghorayeb, K.; Anraku, T.; Firoozabadi, A. Interpretation of the Unusual Fluid Distribution in the Yufutsu Gas-Condensate Field. Soc. Pet. Eng, J. 2003, 8, 114. (53) Hoteit, H.; Santiso, E.; Firoozabadi, A. An Efficient and Robust Algorithm for the Calculation of Gas-Liquid Critical Point of Multicomponent Petroleum Fluids. Fluid Phase Equilib. 2006, 241, 186. (54) Shelton, J. L.; Yarborough, L. Multiple Phase Behavior in Porous Media During CO2 or Rich-Gas Flooding. J. Petrol. Technol. 1977, 29, 1171. (55) Nichita, D. V.; Broseta, D.; Hemptinne, J. C. Multiphase Equilibrium Calculation Using Reduced Variables. Fluid Phase Equlib. 2006, 246, 15. (56) Nichita, V. D.; Broseta, D.; Montel, F. Calculation of Convergence Pressure/Temperature and Stability Test Limit Loci of Mixtures with Cubic Equations of State. Fluid Phase Equlib. 2007, 261, 176.

ReceiVed for reView June 10, 2009 ReVised manuscript receiVed December 16, 2009 Accepted December 28, 2009 IE900948Z

Calculation of Critical Points of Thermodynamic Mixtures with

Recommend Documents