Algebraic Geometric Method for Calculating Phase ... - ACS Publications

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

An algebraic geometric method for calculating phase equilibria from fundamental equations of state Hythem Sidky, Alan C Liddell, Dhagash Mehta, Jonathan D. Hauenstein, and Jonathan K Whitmer Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b02597 • Publication Date (Web): 05 Oct 2016 Downloaded from http://pubs.acs.org on October 6, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35

Su pp

An algebraic geometric method for calculating

or

tin

phase equilibria from fundamental equations of

g

In

state

fo rm

Hythem Sidky,† Alan C. Liddell, Jr,‡ Dhagash Mehta,† Jonathan D.

at

io

Hauenstein,∗,‡ and Jonathan K. Whitmer∗,†

n -F

Department of Chemical and Biomolecular Engineering, University of Notre Dame, Notre

or

Dame, IN 46556, USA., and Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA.

vi

Re

E-mail: [email protected]; [email protected]

ew Abstract

ly

On

Computing saturation properties from highly–accurate Helmholtz equations of state

-N

can be challenging for many reasons. The presence of multiple Maxwell loops often results in incorrect solutions to equations defining fluid phase coexistence. Near the

ot

critical point, the same equations also become ill–conditioned. As a consequence, with-

r fo

out highly accurate initial guesses, it is difficult to avoid the trivial solution. Here, we

b Pu

propose an algorithm applying the technique of Newton homotopy continuation to determine the coexistence curve for all vapor–liquid equilibrium conditions between the

lic

∗

To whom correspondence should be addressed Department of Chemical and Biomolecular Engineering, University of Notre Dame, Notre Dame, IN 46556, USA. ‡ Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA. †

n

io

at

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

1

ACS Paragon Plus Environment


Su

triple point and critical point. Importantly, our algorithm is entirely convergence–

pp

parameter free, does not rely on the use of auxiliary equations, requires no initial guesses, and can be used arbitrarily close to the critical point. It is also fully general-

or

izable to arbitrary equations of state, only requiring that they be locally analytic away

tin

from the critical point. We demonstrate that the method is capable of handling both

g

technical and reference quality fundamental equations of state, is computationally in-

In

expensive, and is useful in both evaluating individual state points and plotting entire

fo

rm

fluid phase envelopes.

at Introduction

n

io -F

The modern fundamental equation of state (FEOS) represents the state-of-the-art in high

or

accuracy thermodynamic property description of fluids. These equations of state are typically explicit in Helmholtz free energy as a function of temperature and density, allowing for

Re

all thermodynamic quantities of interest to be expressed through combinations of partial

vi

derivatives and equality constraints. The ability of FEOS to represent properties within

ew

the error of the most accurate experimental data 1 has driven the development of publicly

On

available software 2,3 supplanting the use of traditional printed tables.

ly

FEOS can be broadly classified into two categories 4 , with reference equations of state

-N

offering superior reliability, extrapolation behavior, and lower uncertainty, particularly near the critical point. The functional forms are often fluid specific and contain a significant num-

ot

ber of terms, with water 5 and carbon dioxide 6 equations containing non-analytic expressions.

r fo

As a consequence, they have been developed for only a handful of additional fluids including

b Pu

nitrogen 7 , argon 8 , methane 9 , ethane 10 , propane 11 , ethylene 12 and sulfur hexaflouride 13 . In contrast, technical equations of state are available for a much wider range of fluids, often

lic

sharing simultaneously optimized forms across entire fluid classes 14–17 . The requirements

at

on experimental data are less stringent with the resulting equations containing exclusively

io

analytic and far fewer terms.

n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 35

2


Page 3 of 35

Su

Due to the explicit Helmholtz energy form, single phase thermodynamic quantities such

pp

as pressure, fugacity, enthalpy, heat capacity and other common properties can be calculated directly at fixed temperature and density by taking the appropriate partial derivatives. On

or

tin

the other hand, in the two–phase region, the equilibrium densities at a fixed temperature must first be identified using an iterative method, after which the desired properties can

g

In

be computed. This can prove to be challenging as FEOS often present multiple solutions

fo

at low temperatures. In addition, the binodal narrows significantly and the Jacobian of

rm

the objective function is poorly conditioned near the critical point. To address these issues,

at

auxiliary equations are fit to saturation data during the development of the initial FEOS and

io

used as a source of high-quality initial guesses. Furthermore, a combination of algorithmic

n

and heuristic approaches assist in convergence and identifying the correct solutions.

-F

As an alternative to the use of successive substitution and classical Newton–Raphson (NR),

or

which are common but require the use of auxiliary equations, Akasaka 18 proposed a damped

Re

NR scheme. In the method, a convergence parameter γ which dampens the Newton step

vi

is adjusted after each iteration. The algorithm’s stability is enhanced and converges suc-

ew

cessfully without the use of auxiliary equations at low to moderate temperatures. However,

On

near the critical point the use of a new starting guess is required. The temperature at which the starting guess is switched depends on the fluid and can range from a few degrees K to

ly

nearly 50 K 18 away from the critical point. Such an approach requires tabulation of limiting

-N

temperatures and is not convenient for maintaining an up–to–date database of fluid equa-

ot

tions. Also, despite being more robust, this method fails within the immediate vicinity of

r fo

the critical point (≈ 0.1 K) due to poor conditioning of the Jacobian 3 . Span 1 utilized the

b Pu

Regula Falsi method, a derivative–free bracketed root finding algorithm, to solve the saturation equations for both pure fluids 1 and multicomponent mixtures 19 . While it overcomes

lic

some difficulties associated with poor conditioning, the rate of convergence is linear and a

at

search interval must be known a priori 20 . This method also proves to be challenging near

3


n

the critical point as the two phase region becomes exceptionally narrow.

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Su

In this work we use numerical algebraic geometry (NAG) to develop a convergence–

pp

parameter free method for determining the saturation states of FEOS. NAG techniques, including various forms of homotopy continuation, have previously been applied to a variety of

or

problems in thermodynamics and statistical mechanics 21–25 including phase stability analysis

tin

of cubic equations of state 26 and critical point determination. 27,28 The algorithm we propose

g

In

here is generalizable to arbitrary EOS models as will be described below, but is specifically

fo

adapted and tested on Helmholtz–explicit FEOS. It uses no auxiliary equations and relies

rm

exclusively on the analytic properties of the system of equations away from the critical point

at

to reliably determine the equilibrium conditions. The use of adaptive multi–precision homo-

io

topy continuation 29 allows the algorithm to robustly obtain solutions arbitrarily close to the

n

critical point where the Jacobian of the objective function becomes ill–conditioned and the

-F

equations non–analytic. Furthermore, utilizing temperature as a natural homotopy param-

or

eter such as this one enables the tracing of the entire binodal and spinodal curves between

Re

the triple and critical points in an efficient manner, using predictor–corrector steps to guide

vi

sampling along the path.

ew On

The Problem Set Up

ly

The modern Helmholtz–explicit FEOS is generally formulated as the sum of the non–

-N

dimensionalized ideal (α0 (τ, δ)) and residual (αr (τ, δ)) Helmholtz energy 1 :

r fo

a(T, ρ) = α(τ, δ) = α0 (τ, δ) + αr (τ, δ) RT

ot (1)

b Pu

where a is the specific Helmholtz free energy, R is the gas constant, T is the temperature, and ρ is the density. The independent variables τ and δ are the inverse reduced temperature

lic

(Tc /T ) and reduced density (ρ/ρc ) respectively. Given the Helmholtz free energy, all ther-

at

modynamic properties in the single phase can be obtained through appropriate derivatives of α. For example, the isochoric heat capacity can be calculated from the reduced Helmholtz 4


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 35

Page 5 of 35

Su

energy as follows, 2 cv 2 ∂ α . = −τ R ∂τ 2 δ

pp or

(2)

tin

The ideal gas contribution α0 is usually defined via correlations of the heat capacity or

g

the Helmholtz free energy directly. The general formulation of the residual term, determined

In

through a complex optimization procedure, can be written as 1

fo

rm

αr (τ, δ) =

at

X

ni τ t i δ d i +

X

X

ni τ ti δ di exp(−γi δ pi )

ni τ ti δ di exp −ηi (δ − εi )2 − βi (τ − γi )2 X + ni δ∆bi exp −ei (δ − 1)2 − fi (τ − 1)2

+

(3)

n

io

where

-F

n 1/(2βi ) o2 ai ∆ = (1 − τ ) + ci (δ − 1)2 + di (δ − 1)2 .

or

(4)

Re

The first two summations are common and found across both technical and reference equa-

vi

tions of state. The Gaussian summation improves accuracy in the critical region and the

ew

final summation is non-analytic and captures the true behavior of the isochoric heat capacity and speed of sound within the immediate vicinity of the critical region. Variables with the

On

subscript i represent parameters that are determined during the fitting process and can be

-N

specific to an individual fluid or a class of substances.

ly

The equilibrium of the liquid phase (′ ) and vapor phase (′′ ) at a given temperature Ts ,

ot

independent of the specific functional form of the equation of state, satisfies the following

r fo

conditions:

p(ρ′ , Ts ) = p(ρ′′ , Ts )

(5)

lic

g(ρ′ , Ts ) = g(ρ′′ , Ts )

b Pu

(6)

io

at where g is the molar Gibbs free energy and p is the unknown saturation pressure. Using the

n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5



Su

simple density dependence of α0 and recognizing that

pp or g

and

p = 1 + δαδr ρRT

(7)

g = 1 + α0 + αr + δαδr RT

(8)

tin In fo

allows for the conditions to be rewritten as

at

rm

δ ′ (1 + δ ′ αr (τ, δ ′ )) = δ ′′ (1 + δ ′′ αr (τ, δ ′′ ))

(9)

-F

and

n

io (10)

δ ′ αδr (τ, δ ′ ) + αr (τ, δ ′ ) + ln δ ′ = δ ′′ αδr (τ, δ ′′ ) + αr (τ, δ ′′ ) + ln δ ′′ .

or

Re

We note that the two simultaneous equations (9) and (10) are implicit in pressure and

vi

have an implicit dependence on the ideal gas free energy α0 through the ln δ terms. With

ew

ρ′ and ρ′′ as the two unknowns, the novel algorithm developed to solve these equations is

On

described below. The obtained densities are then substituted back into (7) to determine pressure, or any other derived thermodynamic identity.

ly -N

We also present an algorithm for situations where the saturation pressure is known and the temperature is unknown. In that case, the condition from equation (10) remains the

ot

same, with two conditions introduced corresponding to the known pressure:

r fo

p = RT δ ′ (1 + δ ′ αr (τ, δ ′ ))

(12)

n

io

at

yielding a system of three equations and three unknowns.

(11)

lic

p = RT δ ′′ (1 + δ ′′ αr (τ, δ ′′ )),

b Pu

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 35

6


Page 7 of 35

Su

Newton homotopy

pp

A core part of the algorithm involves the use of homotopy continuation applied to a Newton

or

homotopy to obtain solutions of a system of nonlinear equations. Mathematically speaking,

tin

a homotopy is a continuous deformation of one continuous function into another. That is,

g

given two functions f, g with domain X and target Y , a homotopy H : X × [0, 1] → Y is a

In

continuous function on X × [0, 1] such that H(x, 0) = f (x) and H(x, 1) = g(x). Numerical

fo

algebraic geometry uses homotopy continuation to solve systems of nonlinear equations. By

rm

casting a target system f as a member of a family of systems containing another system g,

at

one whose solutions are known or readily found, we may continuously deform g into f . The

n

io

solutions of H(x, t) = 0 describe one or more one-real-dimensional curves, called paths which

-F

start at solutions of g = 0 and end at solutions of f = 0.

or

The responsibility of numerical algebraic geometry is tracking these paths, which entails

Re

producing a sequence of numerical approximations to points along these curves, starting at t = 1 and terminating at t = 0. Numerical homotopy continuation uses predictor-corrector

ew

vi

methods to compute these approximations, such as an Euler predictor and a Newton corrector, though higher-order predictors, such as classical Runge–Kutta and its variants, are used

On

in practice. These methods are implemented in the software package Bertini 30,31 , utilized

ly

in this work, which also makes use of adaptive step sizes, adaptive precision 29 , and espe-

-N

cially endgames 32 (see discussion below) to handle the various difficulties that come with

ot

numerically tracking paths, namely regions of ill-conditioning in the vicinity of the critical point.

r fo

Broadly, an endgame is an algorithm which approximates a (possibly singular) solution

b Pu

to H(x, 0) given points along the path. Provided that the solution paths defined by H(x, t) are smooth on (0, 1] and that the system H(x, t) itself is locally analytic for all t ∈ (0, 1],

lic

we may track each solution path from t = 1 to t = 0. Near t = 0, the tracking will employ

at

various endgames to handle singular solutions to H(x, 0) = f (x) = 0. Switching from a path

tracking algorithm to an endgame algorithm is often necessary because, in the ill-conditioned 7


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Su

regions near singularities, pressing on with path tracking can be prohibitively expensive in

pp

terms of precision and step size. We will employ Newton homotopies to force solutions together and to see where they

or

tin

meet. A linear homotopy between f and g has the form

g H(x; t) = (1 − t) · f (x) + t · g(x).

In fo rm

A Newton homotopy 33 is a homotopy in which only the constant terms depend on t. In the context of a linear homotopy, this means that f (x) − g(x) is constant, say v = f (x) − g(x).

at

io

A Newton homotopy has the form

n -F

H(x; t) = f (x) − t · v.

or Suppose that f (x; p) = 0 is a system of equations where x is a variable and p is a

Re

parameter. For a fixed parameter p0 , suppose further that x0 and y0 are two distinct solutions

vi

of f (x; p0 ) = 0. We may wish to determine how to vary the parameter p in order for the

ew

two distinct solution paths starting at x0 and y0 , respectively, to meet at a single solution z ∗

On

corresponding to some new parameter p∗ . Using a Newton homotopy, one may now consider

ly

p as a variable and consider the homotopy: 





-N 

 0  f (x, p)        H(x, y, p; t) :=  f (y, p) − t ·  0  .     x0 − y0 x−y

ot r fo b Pu

The point (x0 , y0 , p0 ) is clearly a solution to the system H(x, y, p; 1) = 0. If the solution path through this point is smooth on (0, 1], then, as t approaches 0, the two distinct solutions

lic

y0 and z0 corresponding to parameter p0 will be forced to merge into a single solution z ∗

at

corresponding to a new parameter p∗ . This observation is a key feature of our methods for

n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 35

computing spinodal and binodal curves. 8


Page 9 of 35

Su

Algorithms

pp

Depending on the information desired from a saturation property calculation, it may be

or

necessary for an underlying software application to implement one or more algorithms. If no

tin

a priori information is known about the system then Algorithm 1 must be used initially to

g

obtain a starting point for the remaining algorithms. Once this starting point is computed,

In

Algorithm 4 may be used to map out the entire spinodal region to the critical point. Algo-

fo

rithm 2 obtains points along the binodals at a fixed temperature. This approach is extremely

rm

useful if the true critical point of the model fluid is unknown. A map of the entire binodal

at

phase diagram can be obtained using Algorithm 3. This algorithm is also necessary if there

n

io

are discontinuities in the spinodals as a function of temperature, as will be described later

-F

on. For applications requiring fixed pressures with unknown saturation temperatures, we provide Algorithm 5.

or Re

The first step in solving the saturation equations is identifying the appropriate spinodals, corresponding to the points where ∂p/∂ρ|T = 0 which are used in subsequent steps. Due to

ew

vi

the lack of constraints on the FEOS between the spinodals, most exhibit undesirable behavior in the unstable two phase region. This often gives rise to multiple Maxwell loops which are

On

problematic for numerical solvers because they produce additional unphysical solutions. For

ly

our purposes, the equations can present additional points satisfying the saturation conditions.

-N

To illustrate this phenomenon, Figure 1 compares the pressure–density isotherms for n-

ot

butane at 300 K using the Peng–Robinson EOS, the original Benedict–Webb–Rubin (BWR) EOS 34 and the technical FEOS by Span and Wagner 15 . The FEOS shows a prominent

r fo

oscillation in contrast to the classical PR–EOS and BWR EOS, which exhibit a single physical

b Pu

Maxwell loop.

Near the critical point, the narrowness of the binodal only complicates matters further. It

lic

becomes challenging to obtain a non-trivial solution without the use of high quality starting

at

guesses. Lemmon and Jacobsen 35 developed a new functional form that yields only one solution for phase equilibrium at a given state, not unlike cubic EOS, through the addition 9


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Su

300

pp or tin g In

100

fo at

P (bar)

rm n

io -F

-100

or Span-Wagner BWR PR

ew

vi

Re

-300

200

400

600

ly

0

On ρ(kg/m3 )

-N

Figure 1: P –ρ isotherm at 300 K for n-butane showing both the metastable and unstable regions using the Span–Wagner technical fundamental equation of state 15 (solid line), Benedict–Webb–Rubin (BWR) equation of state 34 (dashed line) and Peng–Robinson (PR) equation of state (dash-dotted line). Both the BWR and FEOS exhibit additional Maxwell loops which give rise to unphysical solutions.

ot

r fo

b Pu

of certain nonlinear fitting constraints. However, many existing FEOS do not make use of this functional form and its application to reference type FEOS remains limited.

at

lic

In our method, we take advantage of the fact that for many EOS, including FEOS, the outermost stationary points ∂P/∂ρ = 0 always correspond to the true limits of stability 1 .

10


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 35

Page 11 of 35

Su

Algorithm 1: Algorithm for computing spinodals input : Initial temperature T0 , boundary density ρ0 output: Spinodal density ρs that satisfies ∂P =0 ∂ρ T =T0 , ρ=ρs

pp

or

tin

∂P ∂ρ T =T0 , ρ=ρ0

and store result in dP0 ∂P 2. Construct H(ρ; t) := ∂ρ − t · dP0

1. Compute

g

T =T0

In

3. Evolve the homotopy from t = 1 to t = 0 and store result in ρs

fo

rm

These outermost vapor and liquid spinodals can be determined easily and reliably by “rolling

at

downhill” or “climbing uphill” respectively. Any reasonable choice of univariate root finding

io

algorithm is sufficient for solving ∂P/∂ρ|T =T0 = 0 and is initialized by approaching P (ρ; T )

n

from 0 and ∞, which, in practice, is a very large number. We utilize a Newton homotopy

-F

formulated in Algorithm 1 for both 0 and ∞. This step is usually performed at a low

or

temperature, typically near or at the triple point, as there is significant separation between

Re

the liquid–like and vapor-like densities. For EOS such as SAFT and PC–SAFT, the outermost

vi

stationary points may not correspond to the true limits of stability 36 and it is necessary to

ew

employ a procedure other than Algorithm 1 to obtain the desired spinodals. In the case of SAFT–type EOS, a simple algorithm is described in Ref. 36.

On

After computing the spinodals, we present two variants of the algorithm for obtaining

ly

either a solution to the saturation equations at a particular state condition, or alternatively

-N

an entire phase diagram for the fluid. In our first approach, at or near the triple point,

ot

we use the spinodals as starting guesses to compute the vapor and liquid densities. This

r fo

data can be used as the starting point for tracking to a particular temperature using the Newton homotopy defined in Algorithm 2. Rather than tracking to a particular temperature,

b Pu

Algorithm 3 tracks the binodals until the vapor and liquid densities are equal. The advantage

lic

of this algorithm is that an endgame can be used to accurately compute the ending data

at

using data collected along the path potentially far away from the ill-conditioned region near

n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11



Su

Algorithm 2: Algorithm for tracking binodals to a particular temperature input : Start temperature densities ρliq , ρvap satisfying T0 , starting P (ρvap , T0 ) P (ρliq , T0 ) , and target temperature T ∗ . = g(ρvap , T0 ) g(ρliq , T0 ) output: Densities ρ∗liq , ρ∗vap that satisfies P (ρ∗liq , T ∗ ) P (ρ∗vap , T ∗ ) = g(ρ∗vap , T ∗ ) g(ρ∗liq , T ∗ )

pp

or

tin

g

1. Construct a homotopy consistingof Eqns. (9) and (10),with temperature T as a P (ρl , T ) − P (ρv , T )  variable such that H(ρl , ρv , T ; t) := g(ρl , T ) − g(ρv , T )  T − T ∗ − t(T0 − T ∗ ) 2. Evolve the homotopy from t = 1 to t = 0 and store the results in ρ∗l , ρ∗v , Tend . if Tend = T ∗ then return ρ∗liq := ρ∗l , ρ∗vap := ρ∗v else return failure end

In

fo

at

rm

n

io

the critical point.

or

-F

Re

The second approach uses a two-step approach which is summarized in Figure 2. The first

vi

step tracks the entire spinodal curve using the algorithm in Algorithm 4 starting from the

ew

spinodal data computed above at low temperature. In this algorithm, the spinodal densities are tracked as a function of temperature until the critical point is reached, at which point the

On

critical temperature and density is returned. Since the binodal and spinodal curves intersect

ly

at the critical point, the final state of the spinodal homotopy is also a solution to (9) and

-N

(10). At this intersection point, one can use a local tangent cone computation 37 to yield the

ot

first prediction along the binodal curve after which a predictor-corrector tracking employed

r fo

on the Newton homotopy in Algorithm 2 can be used to track down to the low temperature point.

b Pu

While the focus of this work is on temperature–based iterations as they are the most

lic

common, we also provide Algorithm 5 for iterating with respect to pressure. An initial

at

saturation pressure and corresponding densities and temperature must be supplied. These

io

can be readily obtained following the procedure described above for a low temperature

n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 35

12


Page 13 of 35

Su

Algorithm 3: Algorithm for tracking binodals to convergence input : Start temperature densities ρliq , ρvap satisfying T0 , starting P (ρvap , T0 ) P (ρliq , T0 ) = g(ρvap , T0 ) g(ρliq , T0 ) ∂P/∂ρ|T =Tc ρ=ρc = ~0 output: Critical temperature Tc and density ρc satisfying ∂g/∂ρ|T =Tc ρ=ρc

pp

or

tin

g

1. Construct a homotopy consisting of Eqns. (9)and (10), with the density  difference P (ρl , T ) − P (ρv , T ) ρl − ρv as a variable such that H(ρl , ρv , T ; t) :=  g(ρl , T ) − g(ρv , T )  ρl − ρv − t(ρliq − ρvap ) 2. Evolve the homotopy from t = 1 to t = 0 and store the results in ρ∗l , ρ∗v , T ∗ if ρ∗l = ρ∗v then return ρc := ρ∗l , Tc := T ∗ else return failure end

In

fo

at

rm

n

io

-F

solution utilizing Algorithm 4. Using the obtained densities, the saturation pressure can

or

then be calculated directly. Algorithm 5 can then be used for obtaining the saturation

Re

conditions at an alternative pressure. All of the advantages of our algorithms described

vi

for temperature–based iterations are maintained in this case. The only stipulation is that

ew

an initial solution to the equations be provided, as described. A similar homotopy can be

On

constructed if one desires to track the saturation pressure to the critical point. These procedures are not expensive and require only a handful of function evaluations

ly

for most use cases, as detailed below. Both approaches are equally reliable under certain

-N

conditions and discussed in the results. For all homotopies, tracking is performed using the

ot

adaptive step–size Runge–Kutta–Fehlberg (RKF45) predictor–corrector method with multi–

r fo

precision arithmetic and the power series endgame 31,32 . This enables the solver to compute

b Pu

data arbitrarily close to the critical point even with functional forms such as those for CO2 6 and water 5 which contain non–analytic terms since the homotopy only requires the system

lic

to be locally analytic.

at

It is important to note that although the procedures may seem redundant, as the critical

point of a fluid is known beforehand, it is a required step for a truly convergence–parameter 13


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Su

Algorithm 4: Algorithm for tracking spinodals to convergence input : Start temperature T0 , spinodal densities ρliq and ρvap with ρliq 6= ρvap output: Critical temperature Tc , density ρc , satisfying ∂P =0 ∂ρ T =Tc , ρ=ρc

pp

or

tin

1. Construct a homotopy to determine the stationary points of P with respect to ρv,l   ∂P ∂ρ T ρ=ρ1 ∂P  such that H(ρl , ρv , T ; t) := 

g

∂ρ T ρ=ρ2

In

ρl − ρv − t · (ρ∗liq − ρ∗vap ) 2. Evolve the homotopy from t = 1 to t = 0 and store the results in ρ∗l , ρ∗v , T ∗ if ρ∗l = ρ∗v then return ρc := ρ∗l , Tc := T ∗ else return failure end

fo

at

rm

n

io

-F

free method. Technical fundamental equations of state such as those of Span and Wagner 14

or

are not constrained to the supplied critical parameters. Thus, in order to avoid the use of tabulated resulting critical parameters, they can be determined as such, and applicable to

Re

all equations of state. The resulting Tc , ρc of the fluid are then used to initiate the binodal

vi

tracking.

ew

Results

ly

On

The newly developed method is applied to the technical FEOS for nonpolar and weakly polar

-N

fluids 14,15 and reference FEOS for CO2 6 developed by Span and Wagner. The functional form of αr for the former FEOS is

ot

αr (δ, τ ) = n1 δτ 0.25 + n2 δτ 1.125 + n3 δτ 1.5 + n4 δ 2 τ 1.375

r fo

+ n5 δ 3 τ 0.25 + n6 δ 7 τ 0.875 + n7 δ 2 τ 0.625 e−δ 2

+ n8 δ 5 τ 1.75 e−δ + n9 δτ 3.625 e−δ + n10 δ 4 τ 3.625 e−δ

(13)

io

at

lic

3

3

+ n11 δ 3 τ 14.5 e−δ + n12 δ 4 τ 12.0 e−δ ,

2

b Pu

where the parameters ni are fluid specific. This functional form and its polar fluid counterpart 16 have been applied to a broader selection of fluids 17 and are representative of many 14


n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 35

Page 15 of 35

Su

Algorithm 5: Algorithm for tracking binodals to a particular pressure input : Start pressure , temperature  T0 , densities ρliq , ρvap satisfying  P0  P0 P (ρliq , T0 ) , and target pressure P ∗ . P (ρvap , T0 ) =  P0 g(ρvap , T0 ) g(ρliq , T0 ) output: Densities ρ∗liq , ρ∗vap and temperature T ∗ that satisfies     P (ρ∗liq , T ∗ ) P∗ P (ρ∗vap , T ∗ ) =   P∗ ∗ ∗ ∗ ∗ g(ρvap , T ) g(ρliq , T )

pp

or

tin

g

In

fo

1. Construct a homotopy consisting ofEqns. (11) and (12),with pressure p as a P (ρl , T ) − p   P (ρv , T ) − p  variable such that H(ρl , ρv , T, p; t) :=   g(ρl , T ) − g(ρv , T )  p − P ∗ − t(P0 − P ∗ ) 2. Evolve the homotopy from t = 1 to t = 0 and store the results in ρ∗l , ρ∗v , T ∗ , pend . if pend = P ∗ then return ρ∗liq := ρ∗l , ρ∗vap := ρ∗v , T ∗ else return failure end

at

rm

n

io

or

-F

ew

vi

similar 38–42 fluid specific equations.

Re

We first demonstrate the complete mapping of the vapor-liquid phase diagram for n– butane using the proposed algorithm. Figure 2 shows the two–step homotopy process and

On

the resulting diagram. First, the spinodals at the triple point (134.9 K) are determined

ly

using Algorithm 1 by locating the outermost stationary points at fixed temperature. Using

-N

Algorithm 4, the spinodal densities are then tracked to the critical temperature. Once the

ot

critical density and temperature are known, the values are used to initiate tracking of the

r fo

binodal curves following Algorithm 2 back down to the triple point. The result is the entire coexistence curve for the fluid. We emphasize that outside the fluid-specific parameters for

b Pu

αr , no additional information was required.

A good measure of the accuracy for the proposed method is a comparison of the evaluated

lic

critical properties to those reported in the original work 15 . As previously mentioned, this

at

technical FEOS was not constrained to exactly represent the critical point of the fluids, as doing so would require additional terms to fulfill requirements outside the extended critical 15


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Su

Table 1: Critical properties evaluated from the Span–Wagner technical FEOS 15 for nonpolar and weakly polar fluids.

pp

or tin g In fo

Substance Methane Ethane Propane n–Butane n–Pentane n–Hexane n–Heptane n–Octane Argon Oxygen Nitrogen Ethylene Isobutane Cyclohexane SF6

at

rm

n

io

Tc (K) 190.6123 305.5093 369.9388 425.7588 469.6590 507.7945 541.2259 569.5704 150.7914 154.7087 126.2535 282.4544 407.7495 553.6511 318.7239

or

-F

Pc (MPa) 4.6057 4.8871 4.2554 3.8303 3.3688 3.0416 2.7738 2.5066 4.8798 5.0616 3.4045 5.0504 3.6331 3.9836 3.7544

ρc (kg/m3 ) 159.969 196.798 213.325 215.473 235.199 222.825 224.901 227.626 520.530 419.861 307.348 206.281 217.280 262.889 718.886

Re

region 14 . Thus, the critical values used to reduce temperature and density for a given fluid

vi

do not correspond to the same critical values admitted by the FEOS. Table 1 shows the

ew

evaluated critical properties using our method for all fluids in the original work. With the

On

exception of a single fluid, all numbers match those in Ref. 15 identically to numerical

ly

precision. The discrepancy is only in ρc for argon and lies in the second–least significant

-N

digit. There are a number of potential reasons for this, one of which may be the use of higher precision in our calculations; we report our values to an additional significant digit

ot

for reference. Figure 3 shows the reduced saturation curves for all fluids described by the

r fo

technical FEOS.

b Pu

Next, we evaluate the performance and robustness of this method on the Span–Wagner CO2 FEOS 6 that contains a non-analytic term ∆ in (4) which results in a divergence of

lic

the heat capacity at the critical point. It is therefore preferable to adopt Algorithm 3. The

at

reason for this is illustrated in Figure 4. In the neighborhood of the critical point, the

liquid–phase spinodal curve which was tracked from the triple point, exhibits a jump discon16


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 35

Page 17 of 35

Su

tinuity. At ≈ 303.8985 K, the identity of the “outermost” stationary point, and effectively

pp

the Maxwell loop, changes to an interior point which was previously beyond the limit of mechanical stability. The cause of this phenomenon can be understood by looking at ∂P/∂ρ

or

tin

for temperatures near the discontinuity as shown in Figure 5. With increasing temperature, there is an upwards shift in ∂P/∂ρ until the outermost stationary point becomes an inflection

g

In

point. It is doubtful that this discontinuity was intended to capture physical behavior, as

fo

there exists little experimental data this far into the metastable region. More likely is that

rm

the Gaussian damping associated with the non–analytic term for the critical region causes

at

this artifact. This nonetheless highlights additional considerations that may be important

io

when developing new functional forms for fluids.

n

The jump discontinuity observed for CO2 can present difficulties when tracking the spin-

-F

odal curves to the critical point using Algorithm 4. However, using the alternative procedure,

or

Algorithm 3, completely avoids this issue. It is important to emphasize that the previous

Re

scheme is equally valid for systems without jumps in the spinodals, and is quite useful for

vi

plotting the entire saturation curve. Firstly, as with the previous method, the spinodals

ew

at the triple point are identified using Algorithm 1. The spinodal densities are then used

On

as the start system for a homotopy to solve the saturation equations (9) and (10). Once the saturated densities are determined, the binodals are tracked to a pre–specified condition

ly

using Algorithm 2. Figure 6 (a) shows the result of tracking the saturation curve from the

-N

triple point to the critical point. The open circles indicate the points at which the equations

ot

were evaluated. Only 17 evaluations are required to reach a temperature of 303 K, which

r fo

is comparable to the damped Newton method proposed by Akasaka 18 . However, unlike the

b Pu

damped Newton method which fails to converge at temperatures above 303 K, our approach is able to get arbitrarily close to the critical point. This comes at the expense of function

lic

evaluations, which as shown in Figure 6(b), approaches 80 to get within 10−7 of Tc with 10

at

digits of precision on the densities. The reason for this significant increase in function evaluations is evident by looking at the condition number of the Jacobian which is also plotted in

17


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Su

Figure 6(b). There is an exponential increase in the condition number as the critical point

pp

is approached, requiring additional digits of precision to maintain a fixed degree of accuracy. It is unlikely that most scientific applications require such a degree of precision, yet it is still

or

tin

possible using this method. Fewer evaluations are needed to achieve equivalent tolerances for the technical FEOS due to better conditioning near the critical point.

g

In fo

Conclusion

at

rm

We present an algebraic geometric method based on homotopy continuation applied to New-

io

ton homotopies to robustly determine saturation conditions for Helmholtz–explicit funda-

n

mental equations of state. The proposed algorithm does not rely on the use of auxiliary

-F

equations or initial guesses and is entirely parameter–free. Adaptive step–sizes and precision

or

allow the reliable evaluation of saturation properties along the entire saturation curve, and in

Re

particular within the immediate vicinity of the critical point. Two variants of the algorithm

vi

are described in order to deal with potential discontinuities in the spinodal curves, which can be used to map the entire phase diagram.

ew

Though demonstrated for particular equations of state, this approach is fully general

On

and should work with any locally analytic model. Extensions of this method are possible

ly

to multicomponent mixtures where it becomes more challenging to identify the appropriate

-N

densities. Utilizing this method in computer applications can greatly increase the robustness

ot

of density–solvers near the critical point where the trivial solution is difficult to escape, and

r fo

eliminates the need for auxiliary equations.

b Pu

Acknowledgments

lic

All authors greatly acknowledge resources provided by the Notre Dame Center for Research

at

Computing. HS was supported by the National Science Foundation Graduate Fellowship

Program (NSF-GRFP). ACL, DM, and JDH were partially supported by grants from the 18


n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 35

Page 19 of 35

Su

National Science Foundation (ACI-1460032) and the U.S. Army Research Office (W911NF-

pp

15-1-0219) under the Young Investigator Program. JDH also acknowledges support from a Sloan Research Fellowship (BR2014-110 TR14). DM and JKW acknowledge support from

or

tin

startup funds provided by the University of Notre Dame.

g References

In fo

(1) Span, R. Multiparameter Equations of State; Springer Berlin Heidelberg: Berlin, Heidelberg, 2000.

at

rm

io

(2) Lemmon, E.; Huber, M.; McLinden, M. NIST Standard Reference Database 23: Refer-

n

ence Fluid Thermodynamic and Transport Properties-REFPROP. 2013.

-F

(3) Bell, I.; Wronski, J.; Quoilin, S.; Lemort, V. Pure and Pseudo-pure Fluid Thermo-

or

physical Property Evaluation and the Open-Source Thermophysical Property Library

Re

CoolProp. Industrial & Engineering Chemistry Research 2014, 53, 2498–2508.

ew

vi (4) Span, R.; Wagner, W.; Lemmon, E.; Jacobsen, R. Multiparameter equations of state recent trends and future challenges. Fluid Phase Equilibria 2001, 183-184, 1–20.

On

(5) Wagner, W.; Pruß, A. The IAPWS formulation 1995 for the thermodynamic properties

ly

of ordinary water substance for general and scientific use. Journal of Physical and Chemical Reference Data 2002, 31, 387–535.

ot

-N

(6) Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid

r fo

Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa. Journal of Physical and Chemical Reference Data 1996, 25, 1509.

b Pu

(7) Span, R.; Lemmon, E. W.; Jacobsen, R. T.; Wagner, W.; Yokozeki, A. A Reference

lic

Equation of State for the Thermodynamic Properties of Nitrogen for Temperatures

io

at

from 63 . 151 to 1000 K and Pressures to 2200 MPa. Journal of Physical Chemistry

n

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Reference Data 2001, 29 . 19



Su

(8) Tegeler, C.; Span, R.; Wagner, W. A New Equation of State for Argon Covering the

pp

Fluid Region for Temperatures From the Melting Line to 700 K at Pressures up to 1000 MPa. Journal of Physical and Chemical Reference Data 1999, 28, 779.

or

tin

(9) Setzmann, U.; Wagner, W. A New Equation of State and Tables of Thermodynamic

g

Properties for Methane Covering the Range from the Melting Line to 625 K at Pressures

In

up to 100 MPa. Journal of Physical and Chemical Reference Data 1991, 20, 1061–1155.

fo

rm

(10) B¨ ucker, D.; Wagner, W. A reference equation of state for the thermodynamic properties

at

of ethane for temperatures from the melting line to 675 K and pressures up to 900 MPa.

io

Journal of Physical and Chemical Reference Data 2006, 35, 205–266.

n -F

(11) Lemmon, E.; McLinden, M.; Wagner, W. Thermodynamic properties of propane. III. A reference equation of state for temperatures from the melting line to 650 K and pressures

or

up to 1000 MPa. Journal of Chemical and Engineering Data 2009, 54, 3141–3180.

vi

Re

(12) Smukala, J.; Span, R.; Wagner, W. New Equation of State for Ethylene Covering the

ew

Fluid Region for Temperatures From the Melting Line to 450 K at Pressures up to 300 MPa. Journal of Physical and Chemical Reference Data 2000, 29, 1053.

On

(13) Guder, C.; Wagner, W. A reference equation of state for the thermodynamic proper-

ly

ties of sulfur hexafluoride (SF6) for temperatures from the melting line to 625 K and

-N

pressures up to 150 MPa. Journal of Physical and Chemical Reference Data 2009, 38,

ot

33–94.

r fo

(14) Span, R.; Wagner, W. Equations of State for Technical Applications. I. Simultaneously

b Pu

Optimized Functional Forms for Nonpolar and Polar Fluids. International Journal of Thermophysics 2003, 24, 1–39.

at

lic

(15) Span, R.; Wagner, W. Equations of State for Technical Applications. II. Results for

20


n

Nonpolar Fluids. International Journal of Thermophysics 2003, 24, 41–109.

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 35

Page 21 of 35

Su

(16) Span, R.; Wagner, W. Equations of State for Technical Applications. III. Results for

pp

Polar Fluids. International Journal of Thermophysics 2003, 24, 111–162.

or

(17) Lemmon, E.; Span, R. Short Fundamental Equations of State for 20 Industrial Fluids.

tin

Journal of Chemical & Engineering Data 2006, 51, 785–850.

g (18) Akasaka, R. A Reliable and Useful Method to Determine the Saturation State from

In

fo

Helmholtz Energy Equations of State. Journal of Thermal Science and Technology

rm

2008, 3, 442–451.

at

(19) Gernert, J.; J¨ager, A.; Span, R. Calculation of phase equilibria for multi-component

io

mixtures using highly accurate Helmholtz energy equations of state. Fluid Phase Equi-

n

libria 2014, 375, 209–218.

or

-F

(20) Press, W.; Teukolsky, S.; Vetterling, W.; Flannery, B. Numerical Recipes 3rd Edition:

Re

The Art of Scientific Computing, 3rd ed.; Cambridge University Press: New York, NY, USA, 2007.

ew

vi

(21) Bausa, J.; Marquardt, W. Quick and reliable phase stability test in VLLE flash cal-

On

culations by homotopy continuation. Computers & Chemical Engineering 2000, 24, 2447–2456.

ly -N

(22) Malinen, I.; Kangas, J.; Tanskanen, J. A new Newton homotopy based method for the

ot

robust determination of all the stationary points of the tangent plane distance function. Chemical Engineering Science 2012, 84, 266–275.

r fo

(23) Zhang, H. A Review on Global Optimization Methods for Phase Equilibrium Modeling

b Pu

and Calculations. The Open Thermodynamics Journal 2011, 5, 71–92.

lic

(24) Mehta, D.; Chen, T.; Morgan, J. W.; Wales, D. J. Exploring the potential energy

n

io

landscape of the Thomson problem via Newton homotopies. 2015, 142, 194113.

at

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


21



Su

(25) Mehta, D.; Chen, T.; Hauenstein, J. D.; Wales, D. J. Communication: Newton homo-

pp

topies for sampling stationary points of potential energy landscapes. The Journal of chemical physics 2014, 141, 121104.

or

tin

(26) Sun, A.; Seider, W. Homotopy-continuation method for stability analysis in the global

g

minimization of the Gibbs free energy. Fluid Phase Equilibria 1995, 103, 213–249.

In

fo

(27) Wang, M.; Wong, D.; Chen, H.; Yan, W.; Guo, T.-M. Homotopy continuation method

rm

for calculating critical loci of binary mixtures. Chemical Engineering Science 1999, 54, 3873–3883.

at io

(28) Sidky, H.; Mehta, D.; Whitmer, J. Reliable mixture critical point computation using

n

-F

polynomial homotopy continuation. AIChE Journal 2016,

or

(29) Bates, D.; Hauenstein, J.; Sommese, A.; Wampler, C. Adaptive multiprecision path

Re

tracking. SIAM Journal on Numerical Analysis 2008, 46, 722–746.

vi

(30) Bates, D.; Hauenstein, J.; Sommese, A.; Wampler, C. Bertini: Software for Numerical

ew

Algebraic Geometry (2006). Software available at http://bertini. nd. edu

On

(31) Bates, D.; Hauenstein, J.; Sommese, A.; Wampler, C. Numerically solving polynomial systems with Bertini. 2013, 25, xx+352.

ly -N

(32) Sommese, A.; Wampler, C. The Numerical solution of systems of polynomials arising in engineering and science; World Scientific, 2005; Vol. 99.

ot

r fo

(33) Meyer, G. H. On Solving Nonlinear Equations with a One-Parameter Operator Imbedding. SIAM Journal on Numerical Analysis 1968, 5, 739–752.

lic

b Pu

(34) Benedict, M.; Webb, G.; Rubin, L. An Empirical Equation for Thermodynamic Prop-

at

erties of Light Hydrocarbons and Their Mixtures I. Methane, Ethane, Propane and n-Butane. The Journal of Chemical Physics 1940, 8, 334–345.

n

io

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 35

22


Page 23 of 35

Su

(35) Lemmon, E.; Jacobsen, R. A New Functional Form and New Fitting Techniques for

pp

Equations of State with Application to Pentafluoroethane (HFC-125). Journal of Physical and Chemical Reference Data 2005, 34, 69.

or

tin

(36) Privat, R.; Gani, R.; Jaubert, J. N. Are safe results obtained when the PC-SAFT

g

equation of state is applied to ordinary pure chemicals? Fluid Phase Equilibria 2010,

In

295, 76–92.

fo rm

(37) Hao, W.; Hauenstein, J.; Hu, B.; Liu, Y.; Sommese, A.; Zhang, Y.-T. Continuation

at

along bifurcation branches for a tumor model with a necrotic core. J. Sci. Comput. 2012, 53, 395–413.

n

io -F

(38) Lemmon, E.; Huber, M. Thermodynamic Properties of n-Dodecane. Energy & Fuels 2004, 18, 960–967.

or Re

(39) Colonna, P.; Nannan, N.; Guardone, A.; Lemmon, E. Multiparameter equations of state

vi

for selected siloxanes. Fluid Phase Equilibria 2006, 244, 193 – 211.

ew

(40) Ihmels, E.; Lemmon, E. Experimental densities, vapor pressures, and critical point, and

On

a fundamental equation of state for dimethyl ether. Fluid Phase Equilibria 2007, 260, 36 – 48, Jurgen Gmehling Honour Issue.

ly -N

(41) Richter, M.; McLinden, M.; Lemmon, E. Thermodynamic Properties of 2,3,3,3-

ot

Tetrafluoroprop-1-ene (R1234yf): Vapor Pressure and P -ρ-T Measurements and an Equation of State. Journal of Chemical & Engineering Data 2011, 56, 3254–3264.

r fo

(42) Akasaka, R. New Fundamental Equations of State with a Common Functional Form

b Pu

for 2,3,3,3-Tetrafluoropropene (R-1234yf) and trans-1,3,3,3-Tetrafluoropropene (R1234ze(E)). International Journal of Thermophysics 2011, 32, 1125–1147.

n

io

at

lic

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


23



Su pp or tin g In fo at

rm n

io or

-F ew

vi

Re ly

On ot

-N r fo Figure 2: Temperature–density spinodals and binodals for n–butane using the Span–Wagner technical FEOS 14,15 highlighting the novel algorithm where (1) the spinodal is tracked from the triple point to the critical point, then (2) the binodals are tracked back down to the triple point.

n

io

at

lic

b Pu

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 35

24


Page 25 of 35

Su pp or

1

tin g In fo

0.8

at

rm n

io

0.6

Pr

or

-F Re

0.4

ew

vi On

0.2

ly ot

-N

0 0

0.5

1

2

2.5

3

b Pu

1.5 ρr

r fo

Figure 3: Saturation curves for nonpolar and weakly polar fluids 15 obtained using the novel algorithm described in this work.

n

io

at

lic

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


25



Su pp

304.2 ρl ρv

or tin

In 303.8

fo

T (K)

304.0

g at

rm 303.6

n

400

io

450 ρ (kg/m3 )

500

550

or

-F

Figure 4: Spinodal curves near the critical point for CO2 described by the Span–Wagner FEOS 6 showing a jump discontinuity in the “exterior” liquid–phase spinodal.

0.15

ly

0.05

ot

-N

∂P/∂ρ

303.8000 K 303.8985 K 304.0000 K

On

0.1

ew

vi

Re

0

505

510

515

520

530

535

lic

ρ

525

b Pu

-0.05 500

r fo

Figure 5: Density derivative of pressure for CO2 described by the Span–Wagner FEOS 6 . Note the precise temperature at which the outermost stationary point changes identity. The arrows point to the outermost stationary point at various temperatures.

26


n

io

at

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 35

Page 27 of 35

Su pp or tin g In fo at

290

70 14 60

-F

Function evaluations

n

270 260

or

250

16

(b)

io

280

T (K)

80

ρl ρv

Re

240 230

12

50 40

10

30 8 20 6

10

Function evaluations κ(J)

vi

220 300

600 ρ (kg/m3 )

900

1200

0 −102

ew

0

log10 [κ(J)]

(a)

300

rm

−10−1

4 −10−7

−10−4 T − Tc

On

Figure 6: Complete saturation curve for CO2 described by the Span–Wagner FEOS 6 obtained by Newton homotopy tracking of the binodal densities. The points at which the equations are evaluated are shown on the diagram in (a) as open circles. The increase in the number of function evaluations required to achieve a final tolerance of 10−7 at various temperatures approaching the critical point, as well as the condition number κ(J) of the Jacobian are also plotted in (b); see the main text for further discussion.

ly

ot

-N

r fo n

io

at

lic

b Pu

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


27



300

100

P (bar)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

-100 Span-Wagner BWR PR

-300


0

200

400 3

ρ(kg/m )

600

Page 28 of 35

Page 29 of 35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1

0.8

0.6

Pr

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.4

0.2

0


0

0.5

1

1.5 ρr

2

2.5

3

Page 30 of 35

Page 31 of 35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 304.2 50 51 52 53 304.0 54 55 56 57 303.8 58 59 60

T (K)

ρl ρv

303.6


400

450 ρ (kg/m3 )

500

550


303.8000 K 303.8985 K 304.0000 K

∂P/∂ρ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 0.15 50 51 52 0.1 53 54 55 0.05 56 57 58 0 59 60

-0.05 500


505

510

515

520 ρ

525

530

535

Page 32 of 35

Page 33 of 35

T (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

300

ρl ρv

290 280 270 260 250 240 230 220


0

300

600 ρ (kg/m3 )

900

1200


80

16

(b) 70 14 60 12

50 40

10

30

log10 [κ(J)]

Function evaluations

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

8

20 6

10 0 −102

Function evaluations κ(J) ACS Paragon Plus Environment 4

−10−1

−10−4 T − Tc

−10−7

Page 34 of 35

Page 35 of 35Industrial & Engineering Chemistry Research

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25


Algebraic Geometric Method for Calculating Phase ... - ACS Publications

Recommend Documents