On the impact of using volume as an independent variable for the solution of P–T fluid-phase equilibrium with equations of state

The constant pressure–temperature (P–T) flash plays an important role in the modelling of fluid-phase behaviour, and its solution is especially challenging for equations of state in which the volume is expressed as an implicit function of the pressure. We explore the relative merits of solving the P–T flash in two ensembles: mole numbers, pressure and temperature, in which each free-energy evaluation requires the use of a numerical solver; and mole numbers, volume and temperature, in which a direct evaluation of the free-energy is possible. We examine the performance of two algorithms, HELD (Helmholtz free luid-phase equilibria hase stability elmholtz free energy ibbs free energy AFT-VR equation of state energy Lagrangian dual), introduced in Pereira et al. (2012), and GILD (Gibbs free energy Lagrangian dual), introduced here, for the fluid-phase equilibria of 8 mixtures comprising up to 10 components, using two equations of state. While the reliability of both algorithms is comparable, the computational cost of HELD is consistently lower; this difference becomes increasingly pronounced as the number of components is increased. ublis eng–Robinson equation of state © 2014 The Authors. P


Introduction
Equations of state (EOSs) are widely used to represent the physical properties of fluids, including fluid-phase equilibria, in engineering applications (Assael et al., 1996;Kontogeorgis and Folas, 2009). Some EOSs have been shown to reproduce experimental data to a high degree of fidelity, allowing the entire fluid-phase region to be represented in a consistent manner, rendering them powerful modelling tools (Poling et al., 2001;Kontogeorgis and Folas, 2009). While a variety of mathematical forms have been proposed for EOSs, the statistical mechanical formalism used to develop molecular-based equations leads to an EOS written in the canonical ensemble, i.e., with a functional dependence on the variable set consisting of the mole numbers (n), the volume (V), and the temperature (T), and with the Helmholtz free energy as the key thermodynamic function. However, many of the calculations undertaken with EOSs require the variable set n, P (pressure), and T, associated with the constant pressure-temperature (P-T) flash, which is in common use. For example, phase-equilibria data are often reported with P and T as independent variables. Furthermore, the P-T flash is widely used in many process simulation environments. As EOSs become more complex and more computationally demanding to evaluate, the use of appropriate thermodynamic spaces becomes increasingly important to ensure the efficiency and reliability of the computations (Giovanoglou et al., 2009a(Giovanoglou et al., , 2009b).
An important aspect of working with EOSs formulated in terms of the Helmholtz free energy A(n, V, T ) is the route chosen to obtain the isobaric properties. Since the natural ensemble of the Helmholtz free energy is the vector of mole numbers n, volume V, and temperature T, the pressure is obtained as a derivative of A(n, V, T ), i.e., P = −(∂A/∂V ) n,T , which is also a function of the same variable set so that for a given (n, V, T ) it is possible to obtain analytically a unique value for the pressure. However, given P and T, there are often several possible values for the volume and these values (or roots) cannot be obtained analytically when the dependence of the pressure on volume is a polynomial of order greater than three, as in the case for non-cubic equations such as the statistical associating fluid theory (SAFT) (Chapman et al., 1989(Chapman et al., , 1990. Another example of increased computational complexity is the inclusion of an association term, such as in SAFT, in the cubic plus association (CPA) EOS (Kontogeorgis et al., 1996), and in the associated-perturbed-anisotropic chain theory (APACT) Donohue, 1986, 1988;Economou and Donohue, 1990;Elliott et al., 1990). This necessitates the solution of a nonlinear system of equations; although the set of equations has a unique solution (Xu et al., 2002;Kakalis et al., 2006), this leads to an increase in the computational resources required to evaluate the properties of mixtures.
There has been work in the past to assess the computational requirements and efficiency of various EOSs. Mathias  (1986) compared the computational demands of calculating fugacity coefficients and enthalpies for a selection of EOSs available in Aspen Plus (AspenTech, 2014) at the time of their study. They concluded that the calculation of the density (or volume) roots is not the dominant computational cost. Topliss et al. (1988) came to a similar conclusion in their analysis of the same topic in subsequent work. The main aim of our current work is to re-examine this issue, in light of recent developments in the field of EOS modelling, which require the evaluation of more complex algebraic expressions and, consequently, a different partitioning of mathematical operations between the calculation of the density-dependent and density-independent parts of the EOS.
The solution of the P-T flash problem involves a search for the global minimum of the total Gibbs free energy G(n, P 0 , T 0 ), where P 0 and T 0 are the specified pressure and temperature, respectively. The total Gibbs free energy is calculated as a sum over all phases, but the number of phases present at equilibrium is, an unknown quantity a priori. The most stable state of the system may contain a single phase, or multiple phases, and consequently the optimisation problem of directly minimising the system Gibbs free energy over all phases is of nonstandard form: where G T is the total Gibbs free energy of the system, n is a nc × np matrix with the element n i,j representing the number of moles of component i in phase j, np is the number of distinct stable phases, and N G is the maximum possible number of phases present, predetermined by the Gibbs phase rule. Furthermore, i,j is the chemical potential of component i in phase j, n j is the vector of mole numbers in phase j, and n 0 i is the ith element of vector n 0 , i.e., the total number of moles of component i in the mixture.
The solution of problem (1) has been approached in a variety of ways, depending on the particular application and the required balance of computational efficiency and robustness. Here, we focus on techniques that are applicable when no prior knowledge of the phase behaviour of the mixture is available and where more emphasis is therefore placed on robustness. Such an approach is especially appropriate for parameter estimation (Gau et al., 2000) and for molecular design problems (Giovanoglou et al., 2003). In many methods, the unknown number of equilibrium phases is handled through the framework proposed by Michelsen (1982aMichelsen ( , 1982b, in which alternately, np is postulated, a flash calculation is carried out for the fixed number of phases is handled, and then the stability of the postulated solution is tested by minimising the tangent-plane distance function (Baker et al., 1982;Michelsen, 1982a). This process is repeated until a stable solution is obtained. In particular, the problem of identifying instability or metastability has received significant attention. Much of this work has been in the context of applying either deterministic or stochastic global optimisation techniques to minimise the tangent-plane distance function. Some of the many examples of research in this area include work on deterministic algorithms (McDonald and Floudas, 1995;Hua et al., 1996;McKinnon and Mongeau, 1998;Xu et al., 2002;Jalali et al., 2008;Ivanov et al., 2013), and on stochastic schemes (Henderson et al., 2001;Rangaiah, 2001;Teh and Rangaiah, 2002;Bonilla-Petriciolet et al., 2006;Nichita and Gomez, 2009;Srinivas and Rangaiah, 2006;Rahman et al., 2009).
The most popular method for carrying out the flash calculations embedded in this alternating stability test/flash procedure is K-value updating (Rachford and Rice, 1952), which is based on solving the phase-equilibrium equations via successive substitution. While this method is efficient for mixtures exhibiting near-ideal phase behaviour, the algorithm does not necessarily converge when applied to more non-ideal cases, e.g., polymeric systems, as discussed by Heidemann and Michelsen (1995). This failure to converge may be caused by the assumption, made during the solution procedure, that the fugacity coefficients are independent of composition, which becomes increasingly inappropriate as the phase behaviour departs from ideality. An alternative approach is the GFLASH algorithm of Zaydullin et al. (2014), where special care is taken to adopt numerically favorable formulations for phase stability and phase equilibrium, and which has been used to solve large numbers of P-T flash problems in the context of reservoir simulations. Mitsos and Barton (2007) proposed a formulation of the P-T flash problem which is based on obtaining a dual problem from the minimisation of the single-phase Gibbs free energy and the mass-balance constraints. This formulation does not require that the number of equilibrium phases be postulated beforehand and, in addition, has a concave structure that may be exploited to prevent divergence of the calculations, whatever the phase behaviour of the mixture in question. These features make the formulation amenable to the development of algorithms that are generally applicable to any mixture, without the need to introduce assumptions about the underlying phase behaviour.
There are two approaches to tackling the pressure constraint imposed by a P-T flash. The first, and most common, method involves the application of algorithms based directly on the Gibbs free energy at the specified temperature and pressure. This is a function of composition only (G(n; P 0 , T 0 )), but its evaluation at each composition requires the solution of a nonlinear problem to obtain the correct volume root. The second method involves the incorporation of the volume as an explicit variable in the flash problem and consequently working with the function G(n, V ; P 0 , T 0 ) = A(n, V ; T 0 ) + P 0 V so that there is no need to solve for volume roots. Several algorithms have been developed based on this formulation (Nagarajan et al., 1991a(Nagarajan et al., , 1991bXu et al., 2002;Nichita et al., 2006;Pereira et al., 2010Pereira et al., , 2012.
The quantification of the extent to which the need to solve for volume roots when working at fixed pressure affects computational performance is particularly pertinent to P-T flash calculations, since the pressure equation must be solved iteratively many times during such procedures. On the other hand, the addition of volume as an independent variable, in order to avoid this iterative procedure, increases the dimensionality of the Gibbs free energy minimisation problem by one variable, and may therefore lead to an increase in the computational cost.
The focus of our current work is to assess the trade-off, in terms of computational efficiency, between obtaining volume roots through solution of the pressure equation at each Gibbs free energy evaluation and the addition of volume as an explicit variable in the phase-equilibrium calculation, for two duality-based phase equilibrium algorithms. The first, HELD (HElmholtz free energy Lagrangian Dual, Section 2) (Pereira et al., 2012), does not require the solution of the pressure equation; the second, GILD (GIbbs free energy Lagrangian Dual, Section 3), does. In all other respects the two algorithms are similar, and provide an equally high performance in terms of robustness. These algorithms are applied to several case studies to ascertain whether a volume-based formulation offers benefits in terms of either efficiency or robustness. Two EOSs are considered, the statistical associating fluid theory for potentials of variable range (SAFT-VR) (Gil-Villegas et al., 1997;Galindo et al., 1998) and the Peng-Robinson equation (PR) (Peng and Robinson, 1976). The article is organised in the following manner: firstly, the two duality-based algorithms are described, followed by a discussion of the nonlinear solvers employed in GILD to identify the volume roots. The performance of the two algorithms is then compared for a number of systems, exhibiting diverse fluid-phase behaviour, including vapour-liquid equilibrium (VLE), liquid-liquid equilibrium (LLE), and vapour-liquid-liquid equilibrium (VLLE).

The HELD algorithm
The HELD algorithm is described in detail in our previous paper (Pereira et al., 2012) and we highlight only a few aspects here. In the HELD algorithm, a Lagrangian dual problem is solved, derived from the following single-phase formulation: where nc is the number of components in the mixture, A(x, V, T 0 ) is the intensive Helmholtz free energy (in J mol −1 ) expressed in terms of nc − 1 mole fractions, i denotes a specific component in the mixture, x is the vector of nc − 1 component mole fractions x i , such that nc−1 i=1 x i + x nc = 1, X is a bounded set such that X ⊆ [0, 1] nc−1 , V is the molar volume, V and V are, respectively, lower and upper bounds on V (in m 3 mol −1 ), P 0 and T 0 are, respectively, the pressure (in Pa) and temperature (in K) specified for the flash calculation, and x 0 is the vector of nc − 1 component mole fractions, The mole fraction of component nc has been eliminated from the formulation without loss of generality. Problem (2) is a trivial problem with a unique feasible point, but its dual problem possesses valuable properties (Mitsos and Barton, 2007;Pereira et al., 2010). More specifically, its global solution(s) is (are) the stable equilibrium phase(s) of the P-T flash problem at total composition x 0 . The dual problem used takes the form (Pereira et al., 2010(Pereira et al., , 2012 where is the vector of Lagrange multipliers for the mass-balance constraints, Â V ( ) is the dual objective function, and G D is its maximum value. The solution of problem (D x,V ) yields the equilibrium Lagrange multipliers ( * ), and the composition and volume of one HELD is designed specifically for use with EOSs formulated in terms of the Helmholtz free energy. The pressure equation is solved only once in the course of the algorithm, in order to identify a volume root consistent with the feed conditions. Even if several values of the volume can be found at the given P 0 and T 0 , it is not necessary at this stage to identify the most stable root. The identification of an unstable root during this phase does not impair the ability of the algorithm to converge to the correct solution, and where appropriate, to identify the stability of a single-phase solution at the feed composition and to find the corresponding stable volume root at the conditions n 0 , P 0 , T 0 . HELD has been successfully applied to multiphase multicomponent phase-equilibria calculations for numerous mixtures with the PR EOS and with SAFT-type approaches (Artola et al., 2011;Forte et al., 2011;Mac Dowell et al., 2011;Pereira et al., 2012;Rodriguez et al., 2012). Table 1 GILD algorithm for the solution of the P-T fluid-phase equilibrium problem.

Stage I -Stability test and initialisation
Step 1 -Stability test at x 0 Solve stability test up to 10 nc times using tunnelling algorithm. If negative tangent-plane distance found, go to Step 2. If stability test yields stable result, terminate.

Stage II -Identification of candidate stable phases
Step 3 -Solve the outer problem (OP). Set (k) = * and update best upper bound UBD.
Step 4 -Generation of a cutting plane. Solve the inner problem (IP) with fixed (k) from different starting points Add the corresponding variable values x (k) to M.
Step 5 -Search for candidate stable phases.
Update set M * using criteria (PP). If M * contains 2 or more elements then go to Step 7.
Step 6 -Increment iteration counter k = k + 1 and go to Step 3.

Stage III -Acceleration and convergence tests
Step 7 -Minimisation of the Gibbs free energy over all candidate phases.
Solve problem (GMx). If the solution of problem (GMx) is unsuccessful, increment iteration counter k = k + 1 and go to Step 3.

Step 8 -Convergence test
Test consistency of Step 7 Gibbs free energy with best UBD (C1). Test convergence of chemical potentials (C2). If either test fails, increment iteration counter k = k + 1 and go to Step 3.
Step 9 -Check for trace components.
If trace components are present then solve (RT) for their true equilibrium compositions. END GILD;

The GILD algorithm
The GILD algorithm introduced in our current paper is closely related to HELD. The difference lies in the choice of the core thermodynamic function, which in this case is the Gibbs free energy at constant pressure and temperature, which is solely a function of composition. Thus, the problem formulation here is that of Mitsos and Barton (2007), as are the primal and dual problems. Algorithmically, the method is similar to HELD, the main difference being the treatment of the inner problem, in which volume no longer plays an explicit role. For clarity, we briefly outline the main features of the GILD algorithm. In keeping with previous work (Mitsos and Barton, 2007;Pereira et al., 2010Pereira et al., , 2012, the expressions are developed in terms of nc − 1 mole fractions x, rather than n mole numbers n. The algorithmic steps of GILD, as discussed throughout this section, are outlined in Table 1. The primal problem (P) involves the minimisation of a singlephase Gibbs free energy G(x, P 0 , T 0 ), subject to mass-balance constraints, where P 0 is the specified pressure, T 0 is the specified temperature, G(x, P 0 , T 0 ) is the intensive Gibbs free energy at given mole fractions, pressure and temperature, and the other symbols are as defined in Section 2.
A dual of problem (P) can be formed with the nc − 1 massbalance constraints, producing the semi-infinite problem (SIP): Problem (SIP) may be solved with an outer approximation algorithm, in which a linear representation of problem (SIP) is progressively constructed by choosing specific instances of the vector x. This linearised form is known as the outer problem, and yields a guaranteed upper bound UBD on the solution of (SIP), where M is the set of vectors x m , and G P is the value of the Gibbs free energy at the solution of the primal problem, i.e., evaluated at the total composition x 0 , and specified pressure P 0 and temperature T 0 . Although any mole fraction vector would be a valid member of M, it is most efficient to use instances of x m that add the tightest possible constraints to problem (OP), thereby reducing the total number of iterations required for solution of the dual problem. Such constraints are generated by global solution of the inner problem, Problem (IP) is the minimisation of the Lagrangian function with respect to x, for a fixed value of the outer variable vector k . If the global minimum of (IP) is obtained then the solution L (k) is a guaranteed lower bound on the solution of the dual problem. Although the use of global solutions to problem (IP) is likely to result in fewer iterations of the outer approximation algorithm, guaranteed deterministic global optimisation for nonconvex problems is computationally expensive. In GILD we use local minimisation to solve problem (IP). This approach generates useful constraints for the outer problem but does not produce a guaranteed lower bound and therefore GILD may sometimes converge to a metastable solution. Nevertheless, with an appropriate choice of algorithmic parameters, GILD may be shown to have an infinite time guarantee of convergence to the stable solution, as discussed in the context of the HELD algorithm (Pereira et al., 2012). In the remainder of this section, we describe the main steps briefly.

Stage I -stability test and initialisation
The algorithm, outlined in Table 1, begins with a stability test based on the minimisation of the tangent-plane distance function (Michelsen, 1982a) to ascertain whether a phase split occurs at the given conditions (x 0 , P 0 , T 0 ). A tunnelling algorithm (Levy and Montalvo, 1985;Levy and Gomez, 1985) is used, to increase the probability of identifying instability, as compared to the minimisation of the tangent-plane distance function itself. The tunnelling function f T is minimised with respect to the vector of mole fractions x, and is of the form (4) Function f T is very similar to the tunnelling objective function described in Pereira et al. (2012), and the notation from the latter work is retained; d * = d(x * , P 0 , T 0 ) is the value of the best minimum found so far, and x * is the composition vector at this minimum. d(x, P 0 , T 0 ) is the original objective function, which takes the form, where g 0 i is the gradient at the feed conditions (x 0 , P 0 , T 0 ). If a state is identified as being unstable then the GILD algorithm may progress to stage II, beginning with the initialisation of the upper bound UBD and constraint set M.

Stage II -identification of candidate stable phases
If the feed is found not to be stable then the algorithm proceeds to the solution of the dual problem. The outer problem (OP) and the inner problem (IP) are solved alternately in Steps 3 and 4, until the convergence criteria are fulfilled. These convergence criteria are similar to those of HELD (Pereira et al., 2010), except that the derivatives being tested are now those of the Gibbs free energy, rather than the Helmholtz free energy. Additionally, there is no longer a condition on the difference in the volume of the phases, since volume is never considered explicitly in GILD. As in HELD (Pereira et al., 2012), the set M is searched after each major iteration to identify candidate stable phases (Step 5). The mp candidate phases are recorded in a set M * and they must satisfy the following conditions: where UBD is the current lowest upper bound, obtained by solving problem (OP), L (m) is the value of the objective function of problem (IP) corresponding to solution x m and vector (m) , I m is a set of component indices defined by i ∈ I m ⇒ x m i / = x l , where x l is a lower bound on mole fractions, and I m contains the first nc pp = min(nc − 1, 5) components that satisfy this requirement. Thus, nc pp is the number of components to be converged through the solution of the dual problem. Finally, b , and x are user-defined tolerances. Further explanation of these conditions can be found in Pereira et al. (2012). The lower bound on the mole fractions is imposed to avoid numerical difficulties with trace components and is relaxed in a later step in a manner identical to that of HELD (Pereira et al., 2012, Lucia et al., 2000.

Free energy minimisation (Step 7)
We accelerate the approach to a tightly converged solution to the full phase-equilibrium problem once the set M * has at least two members, i.e., two composition vectors. The composition vectors in M * provide a good approximation of the compositions of the stable phases and therefore an excellent initial guess for the direct minimisation of the Gibbs free energy over all phases identified thus far. In GILD, this minimisation takes the form of problem (GM x ), which is in the space of mass numbers (for a total of one mole of mixture): where mp is the number of phases included in the minimisation, G * is the Gibbs free energy at the solution, q is the mp × nc matrix containing the mass numbers in each phase, q i,j is the mass number of component i in phase j, q m i,j is the mass number of component i in phase j in the set M * as produced from stage II, G q (q j , P 0 , T 0 ) is the Gibbs free energy defined in terms of mass numbers, q 0 i is the total mass of component i in the mixture, and MW i is the molar mass of component i in g mol −1 . Formulating this subproblem in mass number rather than mole fraction is favourable because the massbalance constraints are linear, and the problem has better scaling, especially in highly asymmetric mixtures (e.g., polymer-solvent mixtures). If M * does not contain a set of phases that fulfill the mass balance then problem (GM x ) is infeasible (due to the bounds on the mass number variables constraining them to the proximity of the contents of M * ), and the algorithm proceeds to the next iteration via Stage II. If a feasible solution is found, the corresponding mole and phase fractions are stored in a set S * and the algorithm proceeds to Step 8.

Convergence test (Step 8)
The convergence criteria are exactly as described in Pereira et al. (2012) for HELD, and we refer to our earlier work for a discussion of the convergence test. For completeness, the expressions are included here. We undertake two checks in deciding whether to accept the solution of problem (GM x ). Firstly, the following inequality must hold: where g is a convergence tolerance. This ensures that G * is less than, but close to UBD, the 'tentative' (loosely converged) solution to the outer problem, problem (OP). Secondly, the chemical potentials of each component in all phases in S * must agree to within a user-specified tolerance , where i,j denotes the chemical potential of component i in phase j at T 0 and P 0 . If conditions (C1) and (C2) are not satisfied, the iteration counter is incremented and the algorithm returns to Step 3. Otherwise, the contents of S * are considered to be the stable phase simplex, and GILD proceeds to the calculation of trace component compositions in Step 9 if needed.

Refining the compositions of trace components (Step 9)
Any component that is at a lower bound in composition upon convergence of GILD is subject to refinement by problem (RT): where x k,j is the mole fraction of component k in phase j, which is to be determined more precisely, X is the modified search area corresponding to a smaller lower bound (e.g., 10 −50 ), k,j is the chemical potential of the trace component k in phase j, and * k is the equilibrium chemical potential of component k as calculated in the phase in which k is most abundant. The calculation of the chemical potential at the conditions of interest requires the solution of the pressure equation.

Solution of the pressure equation
As mentioned in Section 1, molecular-based EOSs are often formulated in the Helmholtz free energy A(x, V, T ). In this case, pressure P(x, V, T ) is a derivative property, and is obtained from the relationship: In order to obtain the stable volume corresponding to a specified pressure P 0 , one must either find all solutions of the pressure equation P(x, V, T ) = P 0 for given x, T and P 0 , and choose the one corresponding to the lowest Gibbs free energy, or solve a nonconvex optimisation problem. There are a number of different forms that this nonlinear problem (NLP) can take. One option is to solve an unconstrained minimisation of the form where P 0 is the desired pressure. It should be noted that problem (7) may have multiple solutions for a given x, T and P 0 , and indeed, there is no reliable method to predict the number of solutions a priori for higher than cubic EOSs without resorting to guaranteed global optimisation, which can become intractably slow for use in phase-equilibrium calculations. Depending on the conditions, most EOSs have either one or three solutions within the region of physically meaningful volumes for the fluid range, i.e., 0 ≤ Á ≤ 0.74, where Á is the packing fraction (dimensionless volume occupied by the molecules). When three solutions exist, they most often correspond to the liquid, vapour and unstable volumes for the relevant pressure and temperature. The unstable root corresponds to a maximum in the Gibbs free energy, while the liquid and vapour roots are both minima. One of these minima will be local and one global, where the global minimum will correspond to the stable root. In special circumstances, such as at the vapour pressure of a pure component, or an azeotrope, both minima are global, since they represent coexisting phases at the same composition but with different volumes. While many EOSs generally exhibit either one or three volume roots in the physical domain of volume, this is often difficult to prove, and the reliability of many schemes is based on computational experience rather than rigorous mathematical arguments. Some exceptions exist; for example, it was proved (Giovanoglou et al., 2009b) that the augmented van der Waals EOS has three volume roots in the physical region, even though the equation itself is fifth order in volume. In the general case however, the number of solutions present in the physical domain is unknown and there is no mathematical reason to assume that their number would necessarily be limited to three. For example, both Yelash et al. (2005) and Privat et al. (2010) discuss the presence of multiple critical points, and consequently more than three volume roots, in systems represented with the widely used PC-SAFT EOS (Gross and Sadowski, 2001). Privat et al. (2010) further report that the original SAFT (Chapman et al., 1989(Chapman et al., , 1990 and SAFT-VR (Gil-Villegas et al., 1997;Galindo et al., 1998) EOSs do not exhibit such behaviour.
Although the solution of the pressure equation is often required in thermodynamic modelling, the problem is highly nonconvex and care must be taken both when designing algorithms and interpreting results. For example, Lucia et al. (1990) find that the process of obtaining the roots of the Soave-Redlich-Kwong (SRK) EOS (Soave, 1972) through a successive substitution scheme can fall prey to chaotic behaviour. In recent work, Kamath et al. (2010) propose a process-modelling framework for cubic EOS in which constraints are imposed on the volume derivatives of the pressure, to ensure that either a liquid or volume root is obtained as required.
In our current work, two different methods are used to solve the pressure equation in order to investigate the influence of the numerical method chosen for this task on the overall performance of GILD. Both are based on the formulation of problem (7). The first is a quasi-Newton approach (QN) with a line search, based on the method described by Press et al. (2007) and the second is a straightforward Newton-Raphson (NR) method. Both solvers are applicable to any EOS, provided the relevant derivatives are available. The QN method requires only the first volume derivative of the Helmholtz free energy (i.e., the pressure), and the NR method uses both the first and the second volume derivatives of A(x, V, T ). The same convergence criterion is used in both methods: While other methods may be devised to solve the pressure equation efficiently, the use of these two approaches provides useful benchmarks for the purposes of the study undertaken here.

Numerical case studies
In our current work, calculations are undertaken for two equations of state: the Peng-Robinson equation (PR) (Peng and Robinson, 1976) and the statistical associating fluid theory for potentials of variable range (SAFT-VR) (Gil-Villegas et al., 1997;Galindo et al., 1998). The PR EOS has a cubic form in volume, and therefore, for fixed x, P, T the volume roots can be obtained analytically. However, we do not exploit this property here since we use the PR EOS with the aim of providing an assessment of the behaviour of the HELD and GILD algorithms when applied to general EOSs. Clearly, the ability to obtain the pressure roots of an EOS analytically improves the computational efficiency. Indeed, under these circumstances, one would expect GILD to require fewer computational resources, provided that good solutions to the inner problem are added to the constraint set with a similar frequency in both algorithms, i.e., that the rates at which global (or near-global) solutions to the inner problem are located in the two algorithms are comparable.
The HELD and GILD algorithms are implemented in FORTRAN90 and solvers from the Numerical Algorithms Group (NAG) Mark Table 2 Pure component parameters for mixtures represented with the Peng-Robinson (PR) EOS (Peng and Robinson, 1976) where T i c , P i c , and ω i are the critical temperature, critical pressure, and acentric factor of component i respectively. The parameters are as used in Pan and Firoozabadi (1998), Robinson et al. (1978), and Nichita et al. (2002).  Table 3 Binary interaction parameters, k ij , for mixtures represented with the Peng-Robinson (PR) EOS (Peng and Robinson, 1976 (The NAG Library, 2013) are used to solve all sub-problems apart from the solution of the pressure equation (6). The inner problems are solved with the SQP algorithm of Gill et al. (1986aGill et al. ( , 1986b, problem (GM x ) with SNOPT (Gill et al., 2002), and the linear optimisations with a method of Gill and Murray (1978), also described by Gill et al. (1991). Analytical first derivatives are used throughout and analytical second derivatives are used in the Newton-Raphson version of the pressure solver. The bounds on composition are X = [10 −8 , (1 − 10 −8 )] nc−1 , such that the lower bound on mole fraction is x l = 10 −8 . The tolerances are set as: = 0.5, b = 10 −2 , x = 10 −3 , Á = 10 −3 , = 10 −6 , g = 10 −6 , and p = 10 −10 . All runs are performed on an Intel Xeon 2.4 GHz machine running Linux.
We present a number of numerical case studies for each algorithm using both the PR and SAFT-VR EOSs. The parameters for the mixtures modelled with the PR EOS can be found in Tables 2  and 3. The examples modelled with the SAFT-VR SW EOS are the same as those presented in our previous paper (Pereira et al., 2012); the reader is directed to this earlier work for phase diagrams and parameter values for the mixtures in question. We present average statistics of calculations for nine different mixtures, exhibiting varied and often highly non-ideal fluid-phase behaviour. For each mixture, we carry out calculations at 100 different combinations of temperature, pressure and composition, chosen at random and lying between the bounds specified in Tables 4 and 5. Each calculation is repeated 100 times to test reliability, resulting in 10,000 independent calculations for each mixture and algorithmic setup. We consider two repeated calculations to be consistent if the two values of the Gibbs function at convergence, G * , are within 10 −5 of each other. For each mixture, the same 100 state points (temperature, pressure, and composition states) are used for all three algorithmic options (HELD; GILD+NR; GILD+QN) so as to achieve an unbiased comparison. The statistics shown include the average CPU times and number of major iterations required to solve the full phase-equilibrium problem for given P, T and x 0 , as well as the average number of iterations required to solve the pressure equation and the inner problem. The average number of calls to the EOS is also reported. In our implementations of the EOSs, we evaluate the Helmholtz free energy and the corresponding pressure at every call of the EOS, i.e., for given T, V, and x. A more efficient implementation could be developed, thereby reducing the absolute CPU Table 4 Statistics for calculations at randomly selected thermodynamic conditions for systems represented with the PR EOS (Peng and Robinson, 1976). 100 independent repeated calculations are carried out at each of 100 randomly selected unstable compositions, temperatures, and pressures, resulting in 10,000 calculations for each system. The bounds on the random selection of the global composition vector are [10 −4 ,0.999] nc−1 in all cases (where nc i=1 x 0 i = 1). 'Range T 0 ' denotes the bounds on the random selection in temperature, and 'Range P 0 ' the same in pressure. 'IT major' is the number of major iterations (a major iteration involves the solution of both the outer and inner problems), 'CPU' is the number of seconds of single processor computing time, 'Calls to EOS' is the number of times that the equation of state is called, 'IT Vsolve' is the number of iterations required by the volume root solver (only applicable for the GILD statistics), and 'IT inner' is the number of iterations required for the inner problem. All statistics refer to the solution of the entire phase-equilibrium problem and are averaged over the 10,000 calculations. Sd is the corresponding standard deviation. 'GILD + QN' uses a quasi-Newton volume root solver and 'GILD +NR' a Newton-Raphson volume root solver. No failure to obtain the best known solution was observed in any system during the calculations. times reported, but this will not have a significant effect on the relative computational cost, which is the main concern of our current work. In addition to the average statistics reported, we also indicate the standard deviation of each statistic over each set of 10,000 calculations. The mixtures modelled with the PR EOS are: I: a binary mixture of hydrogen-sulphide (H 2 S) and methane (CH 4 ), which exhibits vapour-liquid equilibrium; II: a ternary mixture of methane (CH 4 ), ethane (C 2 H 6 ), and nitrogen (N 2 ), which exhibits vapour-liquid equilibrium; III: a six component synthetic sour gas mixture from Robinson et al. (1978), also studied by Nichita et al. (2002). This mixture contains methane (CH 4 ), ethane (C 2 H 6 ), propane (C 3 H 8 ), hydrogen sulphide (H 2 S), carbon dioxide (CO 2 ), and nitrogen (N 2 ) and exhibits three-phase vapour-liquid-liquid equilibrium (VLLE); IVa: a binary mixture of ethane (C 2 H 6 ) and carbon dioxide (CO 2 ), which exhibits azeotropic vapour-liquid phase equilibrium and is included in order to compare with mixture IVb, which has the same components but is modelled instead with the SAFT-VR EOS.
The parameter values employed to model the interactions between the components in the mixture with the PR EOS are listed in Table 2 and Table 3. The SAFT-VR SW case studies considered include, as well as a binary mixture of ethane and carbon dioxide for comparison with the calculation using PR, four other challenging mixtures that Table 5 Average statistics for calculations at randomly selected thermodynamic conditions for systems represented with the SAFT-VR SW EOS (Gil-Villegas et al., 1997;Galindo et al., 1998). See Table 4 for details of the column headings. No failure to obtain the best known solution was observed in any system during the calculations.  IVb: a binary mixture of ethane (C 2 H 6 ) and carbon dioxide (CO 2 ); V: a ternary mixture of polypropylene glycol (PPG-400, molecular weight 400 g mol −1 ), polyethylene glycol (PEG-600, molecular weight 600 g mol −1 ), and water (H 2 O); VI: a ternary mixture of polypropylene glycol (PPG-400, molecular weight 400 g mol −1 ), polyethylene glycol (PEG-21200, molecular weight 21,200 g mol −1 ), and water (H 2 O); VII: a five component mixture of polyethylene (PE-12000, molecular weight 12,000 g mol −1 ) and light gases: nitrogen (N 2 ), propane (C 3 H 8 ), ethene (C 2 H 4 ), and 1-butene (C 4 H 8 ) and VIII: a ten component mixture of polyethylene (PE-12000, molecular weight 12,000 g mol −1 ), light gases and linear alkanes containing nitrogen (N 2 ), propane (C 3 H 8 ), ethene (C 2 H 4 ), 1-butene (C 4 H 8 ), and the alkanes between n-hexane and ndecane (nC 6 -nC 10 ).

Discussion
The results of the calculations with the Peng-Robinson EOS are shown in Table 4 and those with the SAFT-VR SW EOS in Table 5. The overall trends in computational performance are illustrated in Fig. 1. The statistics over the range of mixtures studied highlight differences between the behaviour of the HELD and GILD algorithms. HELD requires fewer calls to the EOS than GILD in all cases; this is due to the repeated solution of the pressure equation to obtain the volume roots in the GILD algorithm. The difference in the number of EOS calls ranges from a factor of around 5 for some of the two component mixtures, to a factor of 20, for the challenging polymeric mixture VI.
The impact of including the volume as a variable in the inner problem is also characterized by the larger average number of iterations required by HELD to solve this subproblem. The inclusion of volume requires the addition of one extra dimension to the inner problem, and the differences in the number of iterations required by HELD and GILD to solve this problem reflect this. Mixtures I and IV comprise only two components, corresponding to inner problems of one and two dimensions in GILD and HELD, respectively, and show a twofold difference in the iteration count of the inner problem. As the number of components grows, the difference in the dimensionality of the HELD and GILD inner problems remains unchanged at 1 and therefore decreases in relative terms; consequently, the relative difference in the iteration count is reduced. For example, mixture VIII requires around 25 iterations on average for the HELD inner problem, and around 21 iterations for the GILD inner problem.
Both algorithms appear to be reliable at identifying the equilibrium states of the randomly selected unstable conditions, and converge consistently to the best known solutions throughout. A comparison of the two different nonlinear solvers used in GILD to identify volume roots, the quasi and full Newton methods, indicates that in most cases the full Newton solver requires fewer iterations and therefore fewer evaluations of the EOSs.
The robustness of HELD and GILD seems to be comparable in these examples, but HELD requires significantly fewer computational resources to solve the same problems, in terms of CPU time and number of EOS evaluations. This indicates that the extra workload imposed by the larger number of EOS evaluations required by GILD far exceeds any overheads incurred by adding a dimension, namely the volume, to the inner problem in HELD. This is seen to hold for both the PR and SAFT-VR SW EOSs, even though the evaluations with the PR EOS involve significantly fewer floating point operations.
Another interesting feature of the result set is the smaller number of total iterations required by HELD for some of the SAFT-VR SW mixtures, specifically mixtures V, VI, VII and VIII. This is not necessarily expected and perhaps indicates that, in those examples, the solutions to the HELD inner problem are on average better than those of GILD. A common feature of the case studies where this occurs is the presence of a polymeric component. The pressure-volume isotherms for polymer-containing systems can be exceedingly steep functions at high densities, making identification of a stable liquid root challenging. This difference in the number of major iterations is primarily caused by a systematic failure of the nonlinear solver to obtain the stable volume root, at some conditions, since such a failure produces lower quality solutions to the inner problems in which it occurs.

Conclusions
The reliability and relative performance of the Helmholtz free energy Lagrangian dual (HELD) and Gibbs free energy Lagrangian dual (GILD) algorithms for phase-equilibrium calculations with the PR and SAFT-VR EOSs have been investigated. The free energy in the GILD algorithm is obtained via solution of the pressure equation, in order to identify a volume root for a given set of n, P, T conditions. On the other hand, the composition-volume space is used in HELD, so that the Gibbs free energy at the given pressure is obtained at termination without explicitly solving the pressure equation. Both algorithms are found to be of high reliability on the basis of calculations at 100 unique points, each repeated 100 times. The HELD algorithm (Pereira et al., 2012) is found to require less CPU time and fewer EOS evaluations to solve the same problems, a feature which becomes particularly evident as the dimensionality of the mixture or the complexity of the EOSs is increased. In cases where the volume root can be identified in a low-cost manner, for example through the use of the analytical solution of cubic equations for certain EOSs, then it is likely that GILD would outperform HELD. In addition, when dealing with thermodynamic models formulated directly in the n, P, T ensemble, GILD would be the natural choice. However, in other cases, we recommend the use of the HELD approach in the composition-volume space, especially when the complexity of the EOS used and that of the mixture modelled are increased. Finally, although the P-T flash is the prevalent phase-equilibrium formalism in process modelling, the superior performance of HELD raises the question of the effectiveness of the extensive use of the n, P, T ensemble. As higher-than-cubic EOSs framed in terms of the Helmholtz free energy become more widely used, it may be beneficial to adopt formalisms which rely on the more natural n, V, T canonical ensemble.