Development of an Accurate Quintic Equation of State: The Necessity

Subscriber access provided by DEPAUL UNIV

Article

Development of accurate quintic equation of state; On the necessity of the irreducible quadratic polynomial in the attractive term Ju Ho Lee, Sang-Chae Jeon, Jae-Won Lee, Young Hwan Kim, Guen-Il Park, and Jeong Won Kang Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie504118r • Publication Date (Web): 12 Feb 2015 Downloaded from http://pubs.acs.org on February 18, 2015

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 47

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

Development of accurate quintic equation of state; On the necessity of the irreducible quadratic polynomial in the attractive term Ju Ho Lee*a, Sang-Chae Jeona, Jae-Won Leea, Young-Hwan Kima, Guen-Il Parka, Jeong-Won Kangb

a

Nuclear Fuel Cycle Process Development Division, Korea Atomic Energy Research Institute, Yuseong-Gu, Daejeon, 305-353, Republic of Korea b

Department of Chemical and Biological Engineering, Korea University, 145, Anam-ro, Seongbuk-Gu, Seoul, 136-713, Republic of Korea

*Corresponding

author:

Tel.:+82

42

866

6098,

Fax.:+82

[email protected] 1 ACS Paragon Plus Environment

42

868

2605,

E-mail:


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

ABSTRACT

In conventional van der Waals type equation of state (EOS) modeling, the denominator of the attractive contribution is regarded as being a necessarily factorable quadratic polynomial for the closed-form expression of the corresponding Helmholtz free energy. This study evaluates the effect of the opposite case, irreducible quadratic polynomial, on the description of the volumetric behavior of real fluid through a comparison with the critical isotherm data, coexistence curve and PVT isotherm data. In a generalized cubic EOS, an irreducible quadratic polynomial was found to yield an improved description of the flattened region around the critical density while a nonclosed form of Helmholtz free energy is resulted. For an improved volumetric description using irreducible quadratic polynomial and a derivation of the closed-form of Helmholtz free energy, we developed a quintic EOS containing seven parameters, two of which are determined by correlating the critical isotherm data and others by regressing the sub and supercritical properties. In the description of the critical isotherms of nine pure compounds, a comparison with BenedictWebb-Robinson-Soave(BWRS) EOS showed that the present model exhibits slightly larger deviations than BWRS EOS, however, a good agreement was obtained with the reference model of REFPROP 8.0 in the description of the saturated vapor pressure, saturated density and PVT isotherm data over a wide range of temperature.

KEYWORDS: Equation of State; Critical isotherm; Quintic EOS; generalized cubic EOS

2 ACS Paragon Plus Environment

Page 2 of 47

Page 3 of 47

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


1. INTRODUCTION

The design of a chemical process depends heavily on the quality of the required physical properties since they considerably affect the cost evaluation for both the process design and operation.1, 2 As available source for such properties, the equation of state (EOS) provides both phase equilibrium and volumetric properties required for the estimation of the size of various units ranging from pipes to distillation towers. After the beginning of the original van der Waals (vdW) EOS3, a number of EOSs have been developed based on theoretical4-6 or empirical backgrounds7, 8, tested against various experimental data and increasingly used in the chemical industry by being embedded into commercial chemical process simulators. However, studies improving the accuracy of various EOSs are in progress. In developing an accurate EOS, especially for an EOS based on statistical mechanics such as SAFT4 and Lattice-Fluid theory6, an exact reproduction of the critical properties including the critical isotherm was not a major concern owing to a mean-field approximation upon which the models are based. These models showed a good reproducibility of the physical properties at regions far from the critical point, and a satisfactory reproduction of the critical properties was only achieved by combining a complicated crossover theory9; However, in several developments10-13 of SAFT, the exact reproduction of the pure compound critical pressures and temperatures was in fact a major concern and some mean-field based variant models14 were reported to exhibit a reasonable reproduction of the critical behavior. A simultaneous close reproduction of both the critical properties and saturated properties was mainly achieved in EOS based on a virial expansion, categorized into Helmholtz equations.15, 16 Being distinguished from other EOSs, this model employs single or multiple density dependent



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 47

exponential terms enabling both an exact reproduction of the critical constants and a close reproduction of the pure-component critical isotherm. Although they provide a very accurate description of the thermophysical properties over a wide range, however, such accuracy requires many model parameters, increasing the model complexity. Some models were found to yield undesirable multiple Maxwell loops in pure and mixture states and thus the well-developed algorithms implementing a phase equilibrium calculation are also required for the distinction of a stable solution from meta-stable and unstable solutions.15 Van der Waals models, initially started as a cubic EOS, are not conventionally regarded as being applicable to an accurate description of all thermodynamic properties including saturated properties and sub and supercritical volumetric properties. The initial focus on the model improvement was to reproduce exactly the saturated vapor pressure7,

8, 17

, and a reasonable

description of the saturated liquid density was a secondary issue later achieved by modifying the attractive term18 or employing a volume correction strategy19. However, as pointed out by Martin20, the volumetric behavior of real fluids tends to exhibit quadratic behavior implying the limited applicability of all cubic EOSs, and several attempts at overcoming this limit, including quartic21-24 and quintic EOS25, have appeared. The advent of such variants of van der Waals models indicates that the cubic EOS framework should be overcome for the development of accurate van der Waals models, although the framework has provided an excellent background for the development of engineering-based models. One of the important features of van der Waals models lies in the factorable attractive contribution. In cubic26, quartic, and quintic EOS25 modeling, the attractive contribution can be generalized as

Patt = −

a V + m1Vb + m2 b 2

(1)

2


Page 5 of 47

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


where a and b represent the attractive force between molecules and hard core volume, respectively. Various attractive terms27 exist for the combinations of m1 and m2; the condition of m1 = 1 and m2 = 0 corresponds to SRK EOS and that of m1 = 2 and m2 = -1 to PR EOS. Under these conditions, the quadratic polynomial can be factored over the real numbers into a product of linear factors; it becomes v ( v + b ) for SRK EOS and (v + ( 2 + 1)b )( v − ( 2 − 1)b ) for PR EOS. This implies a tacit agreement that the quadratic polynomial in the denominator of eq 1 must be factorable, which is a sufficient condition for the expression of the Helmholtz free energy in a closed-form. This approach28 has also provided an empirical basis for other EOS developments. However, we take issue with the necessity of the factorable quadratic polynomial for an accurate EOS development; This study therefore focuses on two aspects: a desirable EOS framework containing a higher order than the cubic EOS framework, and the effect of factorable and irreducible quadratic polynomials on the performance of the EOS. These issues will be closely examined through a comparison with critical isotherm data. The present model construction follows the procedure repeated in the work of Soave29; First, we present a quintic EOS applicable to the description of a flattened region of critical isotherms and demonstrate the capability of the model by comparing with BWRS EOS29, 30 for critical isotherm data of nine components ranging from Argon to SF6. Then, by assigning appropriate temperature dependence in the model parameters, we demonstrate the performance of the proposed model by comparing with the generated subcritical and supercritical properties reproduced by REFPROP31 for the same systems.



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 47

2. MOEL DEVELOPMENT

2.1 Evaluation of the performance of the generalized cubic EOS in describing critical isotherm data

Although many cubic EOSs have been developed, Kumar and Starling26 presented a generalized cubic EOS containing six model parameters,

(

)

Pb ϕ 1 + d1ϕ + d 2 ϕ2 = RT 1 + d 3ϕ + d 4ϕ 2 + d 5ϕ3

(2)

Now, using the procedure given in Appendix, we relate the model parameters with the critical constants. This requires the assumed form of critical isotherm of eq 2 and the repulsive contribution. The former is expected to have the following form: 3 Pbc F (ϕ / β − 1) = +F RTc 1 + d 3 ϕ + d 4 ϕ 2 + d 5 ϕ 3

(3)

The denominator of the right-side of eq 3 has a cubic density dependence, and thus can be factored into linear and quadratic density-dependent polynomials, 3 Pbc F (ϕ / β − 1) = +F RTc (1 − m0 ϕ ) 1 + m1ϕ + m 2 ϕ 2

(

)

(4)

The value of m0 depends on the repulsive contribution to be employed; 1 for the original vdW EOS, 0.5 for Scott32, and 0.42 for Kim et al.33 The repulsive contribution of the generalized cubic EOS differs to the researches, and this study adopts a repulsive contribution34 assuming a nonlinearity of the free volume

Pbc RT

= rep

ϕ(1 + m3ϕ) 1− ϕ

(5)


Page 7 of 47

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


where subscript rep denotes repulsive contribution, and thus m0 is set to 1. If m3 = 0, eq 5 is reduced to that of the original vdW. As for the denominator of eq 4, conventional cubic EOS modeling favors a quadratic polynomial factorable over the reals since it leads to an analytic expression of the Helmholtz free energy. However, we also consider the quadratic polynomial irreducible over the real numbers, since our concern is to examine the effect of both cases on the volumetric description of real fluids. For the former, eq 4 can be rewritten as 3 Pbc F (ϕ / β − 1) = +F RTc (1 − ϕ ) 1 + m1* ϕ 1 + m 2* ϕ

(

)(

)

(6)

where m1*and m2*are real numbers. From both Constraints introduced in Appendix, we obtain

F=

β 3 + β(m + m2* − 1)

(7)

* 1

(

)(

F (1 − β) = β3 (1 + m3 ) 1 + m1* 1 + m2* 3

)

(8)

For the latter, a negative determinant D = m12 − 4m2 < 0 is obtained and Constraints 1 and 2 lead to

F=

β 3 + β(m1 − 1)

(9)

F (1 − β) = β3 (1 + m3 )(1 + m1 + m2 ) 3

(10)

There are five unknowns, F, β, m1, m2 (or m1* and m2*) and m3, and three constraints Constraints 1, and 2, and the exact reproduction of experimental Zc, exist. Thus two constants should be preferentially given for an evaluation of the remaining model constants and in this study, m2 and m3 are specified. The procedure for obtaining other constants is as follows:



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 47

1) Assume m1 and obtain β by solving eq 10 with F substituted by eq 9. 2) Evaluate the calculated critical compressibility factor by Zccal = F/β 3) If Zccal ≠ Zcexp, adjust the value of m1 and go to step (1) until the equality is made This procedure is also applied to the case of the quadratic polynomial factorable over the reals. Then, through a comparison with the critical isotherm data of methane, we evaluate the performance of the generalized cubic EOS for both cases. This is achieved by optimizing the parameters of m2 and m3 or of m2* and m3* through a minimization of AAD defined by AAD =

1 N

N

∑ 1− i =1

ρ iexp , ρ ical

(11)

where N is the number of data points and ρ is the molar density. In addition to the experimental data35, the critical isotherm generated by REFPROP31 were employed for 0 < ρ ≤ 20 mol/dm3 at 20 intervals since the volumetric data near the critical density are scarce. For the model of eq 4, the optimized parameter values were found to m1 = 1.9068, m2 = 6.3, and m3 = −0.99 yielding AAD = 2.2 %, and for that of eq 6, m1* = 1.0094, m2* = 1.0 and m3* = −0.96 were obtained with AAD = 4.2%. The values of the remaining parameters are listed in Table 1. It is interesting that the optimized values of m1 and m2 were found to produce a negative determinant m12 − 4m2 = −21.564 < 0 although we do not impose the constraint for irreducible quadratic polynomial over the reals in parameter optimization. It seems that as shown in Figure 1, a generalized cubic EOS containing factorable quadratic polynomial describes the overall density region better than that of irreducible quadratic polynomial, especially in a high-density region. However, a close investigation at low and medium density regions reveals that the generalized cubic EOS containing irreducible quadratic polynomial reasonably exhibits a flattened region better than that containing factorable quadratic polynomial.


Page 9 of 47

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


Although the employment of a irreducible quadratic polynomial is found to improve the volumetric behavior around the critical density, this approach yields an undesirable effect; the corresponding derived properties such as the Helmholtz free energy cannot be expressed into a closed-form. The fundamental thermodynamic relation allows residual Helmholtz free energy to be derived from the compressibility factor by ϕ

ϕ

Ares  Z −1 1  Pb =∫ dϕ = ∫  − 1 dϕ RT ϕ ϕ  RT ϕ  0 0

(12)

where subscript res denotes the residual property. The generalized cubic EOS of eq 4 containing an irreducible quadratic polynomial can be alternatively represented as

Pbc ϕ + m3 ϕ 2 ϕ 2 (1 + wϕ) = − a c* RTc 1− ϕ 1 + m1ϕ + m2 ϕ 2

(

)

(

)

(

ϕ + ϕ 2 m1 + m3 − a c* + ϕ 3 m2 + m1 m3 + a c* (1 − w) + ϕ 4 m2 m3 + wa c* = (1 − ϕ) 1 + m1ϕ + m2 ϕ 2

(

)

)

(13)

where ac*=a/b/R/Tc. Then the corresponding residual Helmholtz free energy becomes ϕ 1+ m ϕ  Ares 1 + wϕ 3 dϕ =∫ dϕ − ac* ∫  0 1 + m ϕ + m ϕ2  RT 0 1 − ϕ 1 2  

(14)

It should be noted that for the expression of Ares in a closed form, the numerator in the second integral of the right-side term must be proportional to the first derivative of the denominator, 1 + w ϕ ∝ m1 + 2 m 2 ϕ

(15)

This relation is reduced to w=

2m2 m1

(16)

Since eq 13 is assumed to be a cubic EOS, the quartic term in the numerator should vanish,

w=−

m3 m2 ac*

(17)



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 47

Suppose non-zero m1 and m2 satisfying negative determinant D = m12− 4m2< 0. If the original vdW repulsive contribution is considered (m3 = 0), then w becomes zero by eq 17; however, this value is also inconsistent with the value of w evaluated by eq 16 since eq 16 provides a non-zero w for non-zero m1 and m2. Only eqs 16 and 17 are simultaneously satisfied for the special case of ac* = -m1m3/2. This indicates that until the framework of cubic EOS is conserved, a closed-form of Ares is hard to obtain for EOS containing a irreducible quadratic polynomial

2.2 New quintic EOS and the corresponding attractive term

We have shown that, in Figure 1, a close reproduction of a high-pressure branch and of the flattened region near the critical density cannot be simultaneously achieved under the framework of a cubic EOS. Though the framework enables us to obtain the volume roots conveniently, this should be overcome for a better description of the critical isotherm over the entire region. The scaling hypothesis36 suggests that besides the first and second derivatives being zero, the third and fourth derivatives are also zero, ∂ 3 P / ∂ϕ3

Tc , ϕ=ϕc

= ∂ 4 P / ∂ϕ 4

Tc , ϕ=ϕc

= 0 . This condition

requires at least a quintic EOS and eq 4 can be accordingly modified as follows: Pbc (ϕ / β − 1)5 =F +F RTc (1 − ϕ)(1 + m1ϕ + m 2 ϕ 2 )

(18)

However, as our goal is not to exactly satisfy the scaling hypothesis, we slightly modified the above equation as Pbc (ϕ / β − 1) (1 + c1ϕ + c 2 ϕ 2 ) =F +F RTc (1 − ϕ)(1 + m1ϕ + m2 ϕ 2 ) 3


(19)

Page 11 of 47

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


where c1 and c2 are newly introduced constants. Since the number of model parameters increased, we consider the original vdW repulsive term containing no adjustable constant, m3 = 0, and the application of Constraints 1 and 2 lead to the following relations,

F=

β 3 + β(m1 − c1 − 1)

(20)

F (1 − β ) (1 + c1 + c 2 ) = β 3 (1 + m1 + m 2 ) 3

(21)

Using these relations, eq 19 is rearranged to

(B m − B2 m1 )ϕ + B0 m2 − m22 (F + 1) − B2 2 Pc b B ϕ = − 2 ϕ2 − 1 2 ϕ RTc 1 − ϕ m2 m2 1 + m1ϕ + m2 ϕ 2

(

)

(22)

where

B2 = Fc2 / β 3

( = B + F (1 / β

B1 = B2 + F c1 / β 3 − 3c2 / β 2 B0

1

3

)

− 3c1 / β + 3c2 / β 2

)

(23)

The corresponding residual Helmholtz free energy then becomes ϕ 1 Ares 1 =∫ dϕ − 0 1− ϕ RT m2

1 − m2

∫

ϕ

0

B2 dϕ

 (B1m2 − B2 m1 )ϕ + B0 m2 − m22 (F + 1) − B2  dϕ ∫0  1 + m1ϕ + m2 ϕ2  ϕ

(24)

In spite of irreducible quadratic polynomial, Ares can be expressed in closed form if the numerator in the third integral of the right-side of eq 24 is proportional to the first derivative of the denominator,

B1 m 2 − B 2 m1 : B 0 m 2 − m 22 (F + 1) − B 2 = 2 m 2 : m1

(25)

This is rearranged into the cubic equation of m2 2 (F + 1)m 23 − 2 B 0 m 22 + (2 B 2 + B1 m1 )m 2 − B 2 m12 = 0


(26)


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 47

The residual Helmholtz free energy is then expressed in a closed form

Ares (B m − B2 m1 ) ln 1 + m ϕ + m ϕ2 1   B2ϕ + 1 2 = − ln(1 − ϕ) − 1 2 RT m2  2m2

(

) 

(27)

Equation 26 reduces the number of regressed model parameters by 1 since m2 is obtained by solving this equation. Thus, the parameters to be fitted become c1 and c2, and the remaining parameters can be determined by solving eqs 20, 21, and 26, and Zcexp = Zccal = F/β. As shown in Table 1, the optimized parameter values were c1 = −3.8648 and c2 = 3.7992, resulting in m1 = −1.7162 and m2 = 0.9836 corresponding to the case of an irreducible quadratic polynomial. Figure 2 demonstrates that when compared with eq 4, the critical isotherm reproduced by eq 19 exhibits a better performance near the critical density and a slightly improved performance in the high-pressure region; however, a considerably overestimated pressure is still observed. Another approach should be made for the close reproduction of the critical isotherm over the entire density region including a high-pressure branch. We achieve these requirements by only increasing the degree of freedom by 1. Because the framework of a quintic EOS is to be conserved, the denominator of the right-side of eq 19 is slightly modified as Pbc (ϕ / β − 1) (1 + c1ϕ + c 2 ϕ 2 ) =F +F RTc (1 − ϕ)(1 + m1ϕ + m2 ϕ 2 )(1 + kϕ) 3

(28)

where k is a newly introduced constant. From Constraints 1 and 2, two relations are derived

F=

β 3 + β(m1 + k − c1 − 1)

(29)

F (1 − β ) (1 + c1 + c 2 ) = β 3 (1 + m1 + m 2 )(1 + k ) 3

and accordingly, eq 28 is rearranged into


(30)

Page 13 of 47

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


 Γ  ϕ Pbc Xϕ + Y  = − ϕ 2  0,c + 2  RTc 1 − ϕ  1 + kϕ 1 + m1ϕ + m2 ϕ 

(31)

where

Γ0.c

kB1 − B2 − k 2 B0 + k 3 (Fm1 + m1 + m2 ) = km1 − k 2 − m2

X = Y=

(m1 − k )B2 + B0 m2 k − m2 B1 − k 2 m2 (Fm1 + m1 − m2 ) km1 − k 2 − m2

(

(32)

)

B2 − kB1 + (km1 − m2 )B0 + (F + 1) k 2 m2 − k 2 m12 + m22 − km2 (km1 − m2 ) km1 − k 2 − m2

For the closed-form expression of the Helmholtz free energy, X and Y should be related to

X : Y = 2m2 : m1

(33)

This is rearranged as

(

(

)

) (

k 2 m 2 (F + 1) 2m2 − m12 − m1m 2 + k 2m23 − 2 B1 m2 + B0 m1m 2 + B2 m1

)

+ 2(1 + F )m 23 − 2 B0 m22 + m2 ( B1m1 + 2 B2 ) − m12 B2 = 0

(34)

Then, eq 31 is finally rewritten as follows, Γ (m + 2 m2ϕ )   Γ Pbc ϕ , = − ϕ 2  0,c + 1,c 1 2  RTc 1 − ϕ  1 + kϕ 1 + m1ϕ + m2ϕ 

(35)

where Γ1,c = X/2/m2. It is noteworthy that the constraint of eq 26 is a special case of eq 34 with k = 0, and thus the parameter m2 can be freely assumed since k is introduced to satisfy eq 34. Figure 3 shows that the isotherm reproduced by eq 28, especially at critical density region, is well described, and agrees with the experimental data in a high-pressure region. The resulting attractive term becomes

Γ0,cϕ 2 Γ1,cϕ 2 (m1 + 2m2ϕ ) Patt bc =− − . RTc 1 + kϕ 1 + m1ϕ + m2ϕ 2


(36)


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 47

In Table 1, all model parameters are listed for the different proposed EOS forms, eq 35 (or eq 28), the final version of the proposed quinitc EOS shows the lowest AAD in describing the critical isotherm data of methane.

3. RESULTS AND DISCUSSION

3.1 Optimization of critical parameters

The proposed EOS was tested against the critical P–ρ data of nine components ranging from Ar to SF6. The required model parameters are m1, m2, k, Γ0,c, Γ1,c, and bc; however, three of them are determined by three internal constraints, Zccal = Zcexp, d2P = dP = 0 at φ = β. In this study, m1, m2, and k are regarded as the parameters to be optimized through a minimization of AAD. For propane, similarly with methane, the isotherm data generated by REFPROP besides the experimental data37 were included in the parameter optimization; the data were generated between 0 to 440 kg/m3 at 22 intervals. Other parameters Γ0,c, Γ1,c, and bc are evaluated by the following strategy composed of outer and inner loops; In the outer loop, Γ1,c is adjusted until Zcexp = Zccal is met, and in the inner loop, for the given Γ1,c, m1, m2, and k, β is adjusted until ∂P/∂φ = 0 at φ = φc and T = Tc are made. The remaining parameter Γ0,c is determined in inner loop by rearranging the internal constraint of ∂2P/∂φ2 = 0 at φ = φc and T = Tc,

Γ0,c =

(1 + kβ)3  2

Γ1,c  E 2 2 E E3 4   − 2 + 4 − 2 −   2 3  (1 − β)3 β  G 2 G G G   

(37)

where G = 1 + m1β + m2β2 and E = 2 + m1β. bc is obtained by bc = βvc For an evaluation of the capability of the model, the model was compared with BWRS EOS:


Page 15 of 47

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


(

) (

Z = 1 + e1ρr + e2ρ2r + e3ρ4r + e4ρ2r 1 + e5ρ2r exp − e5ρ2r

)

(38)

where ρr is the reduced density defined by ρr = ρ/ρc. The model parameters are obtained by solving the below equations satisfying both critical constraints and Zccal = Zcexp,

1 + e1 + e2 + e3 + e4 (1 + e5 ) exp(− e5 ) = Z cexp

1 + 2e1 + 3e2 + 5e3 + e4 (3 + 3e5 − 2e52 )exp(− e5 ) = 0

2e1 + 6e2 + 20e3 + e4 (6 + 6e5 − 18e52 + 4e53 )exp(− e5 ) = 0

(39)

where e5 = ψc/(Zcexp)2. In the literature29, ψc and e4 were fixed at 0.06 and 0.5, respectively, and Zc was given by the empirical correlation. However, for an objective comparison, ψc and e4 were refitted to the critical isotherm data of each component, and the experimental values of Zc given in REFPROP were used. In Table 2, the optimized model parameters of both models are listed besides the evaluated AAD for the tested components. Both models were found to have an AAD ranging from 0.61 to 2.36 demonstrating the good applicability of the models toward the description of the critical isotherm over the entire density region. However, as implied in the lower AAD of the BWRS EOS and Figures 4 and 5, BWRS EOS reproduces the overall isotherm data better than the proposed model;the flattened region of the experimental data more closely agrees with that reproduced by BWRS EOS than that by eq 35. A reduction of the number of the model parameters is desirable in EOS modeling. From the optimized values of the model parameters listed in Table 2, we found that without the deterioration of the accuracy, k can be related to m1 and m2 by k = −0.8226 + 0.1722 m1 + 0.1753m2

(40)

In Table 3, the newly fitted m1 and m2 are listed besides other constants determined by solving the critical conditions and Zcexp = Zccal.



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 47

The good performance of BWRS EOS was explained by investigating the characteristic behavior of the exponential terms given in eq 39, only active at the intermediate density region.38 Our model can also reproduce a flattened region close to that of the experimental data if the quadratic polynomial in the numerator of eq 28 can be approximated as follows: 1 + c1ϕ + c 2 ϕ 2 ≅ (ϕ / β − 1)

2

(41)

However, as Figures 4 and 5 illustrate, it seems impossible to satisfy this approximation with the current EOS formulation.

3.2 Temperature dependency in model parameters

Once the critical parameters are determined, the appropriate temperature dependence in model parameters allows the description of the thermodynamic properties at sub and supercritical region. Among the model parameters, the temperature dependence is only considered in Γ0, Γ1, and b since m1, m2, and k are directly coupled with the reduced density. From the critical isotherm equation represented by eq 35, the generalized form at the given temperature then becomes Γ ϕ2 Pb Γ (m + 2m2ϕ )ϕ 2 ϕ = − 0 − 1 1 , RT 1 − ϕ Tr 1 + kϕ Tr 1 + m1ϕ + m2ϕ 2

(42)

where Γ0,c = Γ0(Tc), Γ1,c = Γ1(Tc) and bc = b(Tc). Here Tr means the reduced temperature defined by Tr = T/Tc. The temperature dependence for each parameter can be determined through many trials by correlating various thermodynamic properties, however, we found that the behavior of each contribution constituting eq 36 should be preferentially understood. Let us consider the attractive contribution of methane calculated with the critical parameters listed in Table 3. In Figure 6, the 16 ACS Paragon Plus Environment

Page 17 of 47

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


red dashed line exhibiting a negative pressure represents the first term on the right side of eq 36, the orange dashed line the second term and the red solid line the total attractive term evaluated by summing both terms. Noteworthy in this figure is the characteristic behavior of the second term; it roughly exhibits a positive pressure contribution except at the narrow region φ < -2m1/m2 in which a negative pressure contribution is found. The opposite sign of the first and second term causes no problematic behavior at a critical temperature over the entire density region; the magnitude of the last term is much smaller than that of the second term, resulting in negative pressure for the attractive term. However, for certain temperature dependence in Γ0 and Γ1, an undesirable phenomenon may arise at a low temperature region. Here we consider the isotherms of methane reproduced at Tr = 0.3 with the same parameters used in the construction of the previous plot. The temperature, although below the triple point, is chosen only for the purpose of illustration. We suppose two cases of temperature dependent Γ1 :A) Γ1 = Γ1,c / Tr0.5 and B) Γ1 = Γ1,c / Tr2.5. In both cases, Γ0 is assumed to have a temperature dependence of Γ0 = Γ0,c / Tr. Figure 7a represents the isotherms and constituting contributions containing the temperature dependence of case A and Figure 7b to case B. The black solid line indicates the total pressure and the blue line the pressure contribution of the repulsive contribution. In Figure 7a, the overall pressure isotherm demonstrates no particular behavior, however, in Figure 7b, another Maxwell loop is shown in the isotherm. This undesirable behavior originates from the competitive effect between the last term and the second term of eq 36 which show opposite sign; the magnitude of the orange dashed line is similar to that of the red dashed line yielding a weakened attractive contribution compared with that of Figure 7a. In the center of this behavior is the considerably increased Γ1, the coefficient of the last term of eq 36; Γ0 = 11.41 >> Γ1 = 0.625 for case A and Γ0 = 11.41 > Γ1



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 47

= 6.94 for case B. Therefore, to avoid additional Maxwell loops at a low temperature region, among the terms constituting the attractive contribution, the model parameters coupled to the term exhibiting a positive sign must have a weak temperature dependence compared with those exhibiting a negative sign and we consider the following empirical forms for Γ0 and Γ1;   κ Γ0 = Γ0 ,c 1 + κ 1 − κ 2 − 1 + κ 2 Tr  Tr  

(43)

  ε Γ1 = Γ1,c 1 + ε 1 − ε 2 − 1 + ε 2 Tr    Tr  

(44)

Note that in Γ1, compared with 1/ Tr in Γ0, the weak temperature dependence represented by 1/Tr0.5 is employed. For a hard core volume, we assumed the following empirical temperature dependent form;  0 .5  b = bc 1 + δ exp( − )  / (1 + δ exp( −0.5) ) Tr  

(45)

where δ is component-specific parameter and bc = b(Tc) is satisfied. Figure 8 shows the predicted isotherms of CO2 at three temperatures, Tc/2, Tc/5 and Tc/10. In all isotherms, only one Maxell loop representing a vapor-liquid transition is shown, supporting the validity of the temperature dependency introduced in eqs 43 and 44. This feature also allows a fast and robust evaluation of the density root at the given pressure because there is no need for checking the additional Maxwell loop. Although an algorithm52 is available for solving a quintic EOS, the numerical solving for the density root53 is adopted for the present model.

3.3 Extension to subcritical and supercritical region.


Page 19 of 47

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


The determination of five constants, ε1, ε2, κ1, κ2, and δ, requires accurately measured thermodynamic data ranging from saturated properties to PVT data. For this purpose, rather than using limited data sources, we used NIST ThermoData Engine54, 55 (TDE) designed to perform critical data evaluation dynamically “on demand” cased on available experimental data and their uncertainties. In parameter optimization, excluded are the data points in which TDE shows one of following statuses, a large deviation, out of range, invalid value, and predicted values. All constants are simultaneously determined by minimizing the objective function defined as a sum of AAD of the experimental data as follows; cal cal ρsat, ρsat, 100 N1 Pi sat,cal 100 N2 100 N3 l,i v,i OBJF = 1 − sat,exp + ∑ 1 − Psat,exp + N ∑ ∑ 1 − ρsat,exp N1 i=1 ρl,i N 3 i=1 i 2 i =1 v,i

ρcal ρcal 100 N4 100 N5 v, j l, j + 1 − exp + 1 − exp ∑ ∑ N 4 j=1 ρ v, j N 5 j=1 ρl, j

(46)

where N1 , N2, and N3 are the number of experimental saturated vapor pressure, liquid density and vapor density respectively. N4 indicates the total number of experimental PVT data of a single vapor phase at a subcritical region and of a supercritical phase, however, N5 is the number of experimental PVT data of only a single liquid phase at the subcritical region. Superscript sat denotes the saturated properties and subscripts v and l denotes the vapor and liquid phases respectively. The optimized constants of eqs 43-45 for the tested molecules are listed in Table 3 . In demonstrating the performance of the proposed model, instead of the direct comparison with experimental data, we choose to compare the model with an existing accurate model, the reference model of REFPROP, because many data points of different literatures are employed in our model parameter optimization. For a comprehensive evaluation, besides the saturated properties including vapor pressure, liquid density and vapor density, PVT isotherms at several subcritical and supercritical regions are also reproduced by REFPROP, which are shown as the 19 ACS Paragon Plus Environment


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 47

empty circles in Figure 9. A good agreement between both models is found for methane and for other components as revealed in Supporting Information. As revealed in Figure S2, slightly disagreement is shown in the vapor pressure of propane below 140K and the PVT behavior of the liquid phase at high density region. The former can be improved if more temperature dependent terms are considered in the attractive parameter given by eq 43, however, the latter is principally associated with the density dependence of vdw repulsive term. Figures S1 –S9 shows that the liquid density predicted by our model has a steeper slope than the density reproduced by REFPROP, implying a strong density dependency of the original vdw repulsive term as the density is increased. A similar behavior has been also reported in a comparison of vdw repulsive term with the molecular simulation data.27 Consequently, this behavior can also be improved if a weaker repulsive term such as that of lattice fluid is considered. All of these results indicate the applicability of the present approach toward the description of the thermodynamic properties of a pure compound.

4. CONCLUSION

This study is distinguished from other EOS modelings in quantitatively evaluating the effect of an irreducible quadratic polynomial in an attractive contribution on the reproduction of a critical isotherm of the pure components. Although a developed quintic EOS containing the proposed term is slightly less accurate than BWRS EOS in correlating the critical isotherm data, the model yields no additional Maxwell loops in the isotherms and a comparison with the reference model of NIST shows the good applicability of the present model in describing subcritical and supercritical properties over a wide range of temperature.


Page 21 of 47

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


One of the advantages of the present study lies in the generation of additional Maxwell loops in the pressure isotherm at low temperature. Since only single loop exists, the density-root finding is easily implemented and thus the computational load can be significantly reduced compared with the models exhibiting multiple Maxwell loops.56-58 In obtaining thermodynamic properties, an accurate EOS development generally prefer to derive a model by proposing a model for the Helmholtz free energy, and thus other thermodynamic properties such as the pressure may be easily derived through differentiation. In this aspect, the approach followed by the present study, the derivation of the Helmholtz free energy by proposing a model for EOS and integrating the model, may not be straightforward as the conventional approach; however, for an EOS containing a few terms, it allows a direct evaluation of the capability of a given model in describing a critical isotherm through an algebraic manipulation. The role of the irreducible quadratic polynomial is similar to that of the exponential terms of BWRS EOS in improving the experimental volumetric behavior around the critical density. Although its performance is slightly lower than that of the exponential term, the newly proposed attractive term has the advantage of being readily combinable with other repulsive terms of statistical-based models6, 59 because the attractive contribution of these models, the original van der Waals attractive term, can be replaced with the newly derived attractive contribution by eq 36. A further study is in progress to evaluate the applicability of the proposed attractive term into other repulsive terms.

5. ACKNOWLEDGMENT



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 47

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST)(No. 2012M2A8A5025696)

NOMENCLUTURE a

attractive energy parameter

ac*

reduced attractive energy parameter at critical temperature

A

Helmholtz free energy

b

hard core volume

Bi

constants defined by eq 23

ci

constants of eq 19

di

the model parameters of generalized cubic EOS given by eq 2

D

determinant

ei

the parameters of BWRS EOS

F

the reduced critical pressure at critical density given by F = Pcb/R/Tc

k

the parameter of eq 19

m1 , m2

the coefficients of quadratic polynomial

m1*, m2*

the coefficients of quadratic polynomial factorable over the reals

m3

constants in eq 5

N

the number of data point

P

pressure

R

gas constant

T

absolute temperature


Page 23 of 47

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


v

molar volume of the fluid

X,Y

variables given by eq 32

Z

compressibility factor

Subscript att

attractive contribution

c

critical point

l

liquid phase

r

reduced property

rep

repulsive contribution

res

residual contribution

v

vapor phase

Superscripts cal

calculated property

exp

experimental property

sat

saturated property

Greeks φ

reduced volume fraction

β

reduced critical volume fraction

ρ

molar density

ρr

reduced molar density defined by ρr =ρ/ρc



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, Γ1 attractive energy parameters for eq 36 δ

constant in eq 45

ε1, ε2 constants in eq 44 κ1, κ2 constants in eq 43


Page 24 of 47

Page 25 of 47

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


REFERENCES

1.

Cox, K. R., Physical Property needs in Industry. Fluid Phase Equilib. 1993, 82, (0), 1526.

2.

Gupta, S.; Olson, J. D., Industrial Needs in Physical Properties. Ind. Eng. Chem. Res.

2003, 42, (25), 6359-6374. 3.

Waals, J. D. v. d. Over de continuiteit van den gas- en vloeistoftoestand. A.W. Sijthoff, Leiden, 1873.

4.

Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M., New reference equation of state for associating liquids. Ind. Eng. Chem. Res. 1990, 29, (8), 1709-1721.

5.

Vera, J. H.; Prausnitz, J. M., Generalized van der Waals theory for dense fluids. Chem. Eng. J. 1972, 3, (0), 1-13.

6.

Sanchez, I. C.; Lacombe, R. H., An elementary molecular theory of classical fluids. Pure fluids. J. Phys. Chem. 1976, 80, (21), 2352-2362.

7.

Soave, G., Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 1972, 27, (6), 1197-1203.

8.

Peng, D.-Y.; Robinson, D. B., A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, (1), 59-64.

9.

Kiselev, S. B.; Ely, J. F., Crossover SAFT Equation of State: Application for Normal Alkanes. Ind. Eng. Chem. Res. 1999, 38, (12), 4993-5004.

10.

Chen, J.; Mi, J.-g., Equation of state extended from SAFT with improved results for nonpolar fluids across the critical point. Fluid Phase Equilib. 2001, 186, (1–2), 165-184.

11.

Cismondi, M.; Brignole, E. A.; Mollerup, J., Rescaling of three-parameter equations of state: PC-SAFT and SPHCT. Fluid Phase Equilib. 2005, 234, (1??), 108-121. 25 ACS Paragon Plus Environment


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

12.

Page 26 of 47

Polishuk, I., Hybridizing SAFT and Cubic EOS: What Can Be Achieved? Ind. Eng. Chem. Res. 50, (7), 4183-4198.

13.

Polishuk, I., Standardized Critical Point-Based Numerical Solution of Statistical Association Fluid Theory Parameters: The Perturbed Chain-Statistical Association Fluid Theory Equation of State Revisited. Ind. Eng. Chem. Res. 53, (36), 14127-14141.

14.

Beckman, E. J.; Porter, R. S.; Koningsveld, R., Mean-field lattice equations of state. 2. The role of the segmental contact surface area in the unification of recent lattice models. J. Phys. Chem. 1987, 91, (25), 6429-6439.

15.

Span, R., Multiparameter Equations of State: An Accurate Source of Thermodynamic Property Data. Springer: 2000.

16.

Kunz, O.; Wagner, W., The GERG-2008 Wide-Range Equation of State for Natural Gases and Other Mixtures: An Expansion of GERG-2004. J. Chem. Eng. Data 2012, 57, (11), 3032-3091.

17.

Twu, C. H.; Bluck, D.; Cunningham, J. R.; Coon, J. E., A cubic equation of state with a new alpha function and a new mixing rule. Fluid Phase Equilib. 1991, 69, (0), 33-50.

18.

Fuller, G. G., A Modified Redlich-Kwong-Soave Equation of State Capable of Representing the Liquid State. Ind. Eng. Chem. Fundam. 1976, 15, (4), 254-257.

19.

Péneloux, A.; Rauzy, E.; Fréze, R., A consistent correction for Redlich-Kwong-Soave volumes. Fluid Phase Equilib. 1982, 8, (1), 7-23.

20.

Martin, J. J., Cubic equations of state - Which? Ind. Eng. Chem. Fundam. 1979, 18, (2), 81-97.

21.

Soave, G., A non-cubic equation of state for PVT and phase equilibrium calculations (pure compounds). Fluid Phase Equilib. 1990, 56, (0), 39-57.


Page 27 of 47

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


22.

Kubic Jr, W. L., A quartic hard chain equation of state for normal fluids. Fluid Phase Equilib. 1986, 31, (1), 35-56.

23.

Shah, V. M.; Bienkowski, P. R.; Cochran, H. D., Generalized quartic equation of state for pure nonpolar fluids. AlChE J. 1994, 40, (1), 152-159.

24.

Zhi, Y.; Meiren, S.; Jun, S.; Lee, H., A new quartic equation of state. Fluid Phase Equilib.

2001, 187–188, (0), 275-298. 25.

Koziol, A., Quintic equation of state for pure substances in sub- and supercritical range. Fluid Phase Equilib. 2008, 263, (1), 18-25.

26.

Kumar, K. H.; Starling, K. E., The most general density-cubic equation of state: application to pure nonpolar fluids. Ind. Eng. Chem. Fundam. 1982, 21, (3), 255-262.

27.

Wei, Y. S.; Sadus, R. J., Equations of state for the calculation of fluid-phase equilibria. AlChE J. 2000, 46, (1), 169-196.

28.

Yun, Z.; Qin, J.; Shi, M.; Shi, J., A New General Equation of State. Chin. J. Chem. Eng.

1999, 7, (2), 161-165. 29.

Soave, G. S., A Noncubic Equation of State for the Treatment of Hydrocarbon Fluids at Reservoir Conditions. Ind. Eng. Chem. Res. 1995, 34, (11), 3981-3994.

30.

Soave, G. S., An effective modification of the Benedict–Webb–Rubin equation of state. Fluid Phase Equilib. 1999, 164, (2), 157-172.

31.

Lemmon, E.; Huber, M. L.; McLinden, M. O. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties - REFPROP, 8; 2007.

32.

Scott, R. L., In Physical Chemistry, An Advanced Treatise. Academic Press: 1971; Vol. Vol. 8A.



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

33.

Page 28 of 47

Kim, H.; Lin, H. M.; Chao, K. C., Cubic chain-of-rotators equation of state. Ind. Eng. Chem. Fundam. 1986, 25, (1), 75-84.

34.

Hajipour, S.; Edalat, M., A New Hard Sphere Cubic Equation of State for Predicting Fluids’ Properties and Vapor-Liquid Phase Equilibrium Calculations. J. Phs. Eqil. and Diff. 2008, 29, (4), 322-332.

35.

Händel, G.; Kleinrahm, R.; Wagner, W., Measurements of the (pressure, density, temperature) relation of methane in the homogeneous gas and liquid regions in the temperature range from 100 K to 260 K and at pressures up to 8 MPa. J. Chem. Thermodyn. 1992, 24, (7), 685-695.

36.

Eubank, P. T.; Hall, K. R., New Identities for Critical-Point Constraints with Application to Cubic Equations of State. Ind. Eng. Chem. Res. 2004, 43, (15), 4446-4451.

37.

Glos, S.; Kleinrahm, R.; Wagner, W., Measurement of the (p, ρ,T) relation of propane, propylene, n-butane, and isobutane in the temperature range from (95 to 340) K at pressures up to 12 MPa using an accurate two-sinker densimeter. J. Chem. Thermodyn.

2004, 36, (12), 1037-1059. 38.

Kedge, C. J.; Trebble, M. A., Improvements to a new equation of state for pure components. Fluid Phase Equilib. 2004, 215, (1), 91-96.

39.

Funke, M.; Kleinrahm, R.; Wagner, W., Measurement and correlation of the (p, ρ, T) relation of ethane I. The homogeneous gas and liquid regions in the temperature range from 95 K to 340 K at pressures up to 12 MPa. J. Chem. Thermodyn. 2002, 34, (12), 2001-2015.

40.

Douslin, D. R.; Harrison, R. H., Pressure, volume, temperature relations of ethylene. J. Chem. Thermodyn. 1976, 8, (4), 301-330.


Page 29 of 47

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


41.

Nowak, P.; Kleinrahm, R.; Wagner, W., Measurement and correlation of the (p, ρ,T) relation of ethylene I. The homogeneous gaseous and liquid regions in the temperature range from 105 K to 340 K at pressures up to 12 MPa. J. Chem. Thermodyn. 1996, 28, (12), 1423-1439.

42.

Gilgen, R.; Kleinrahm, R.; Wagner, W., Measurement and correlation of the (pressure, density, temperature) relation of argon I. The homogeneous gas and liquid regions in the temperature range from 90 K to 340 K at pressures up to 12 MPa. J. Chem. Thermodyn.

1994, 26, (4), 383-398. 43.

Nowak, P.; Kleinrahm, R.; Wagner, W., Measurement and correlation of the (p,ρ,T) relation of nitrogen. I. The homogeneous gas and liquid regions in the temperature range from 66 K to 340 K at pressures up to 12 MPa. J. Chem. Thermodyn. 1997, 29, (10), 1137-1156.

44.

Duschek, W.; Kleinrahm, R.; Wagner, W., Measurement and correlation of the (pressure, density, temperature) relation of carbon dioxide I. The homogeneous gas and liquid regions in the temperature range from 217 K to 340 K at pressures up to 9 MPa. J. Chem. Thermodyn. 1990, 22, (9), 827-840.

45.

Michels, A.; Michels, C., Isotherms of CO2 between 0 degrees and 150 degrees and Pressures from 16 to 250 Atm (Amagat Densities 18-206). Proc. R. Soc. London, Ser. A

1935, 153, (878), 201-214. 46.

Michels, A.; Blaisse, B.; Michels, C., The Isotherms of CO2 in the Neighbourhood of the Critical Point and Round the Coexistence Line. Proc. R. Soc. London, Ser. A 1937, 160, (902), 358-375.



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

47.

Page 30 of 47

Kang, T. L. Thermodynamic Properties of Sulfur Dioxide. Ph. D. dissertation in chemical engineering, University of Texas in Austin, 1960.

48.

Kang, T. L.; Hirth, L. J.; Kobe, K. A.; McKetta, J. J., Pressure-Volume-Temperature Properties of Sulfur Dioxide. J. Chem. Eng. Data 1961, 6, (2), 220-226.

49.

Hellwig, L. R. Pressure-Volume-Temperature Properties of Sulfur Dioxide. Ph. D. dissertation in chemical engineering, University of Texas in Austin, 1955.

50.

Funke, M.; Kleinrahm, R.; Wagner, W., Measurement and correlation of the (p, ρ, T) relation of sulphur hexafluoride (SF 6). I. The homogeneous gas and liquid region in the temperature range from 225 K to 340 K at pressures up to 12 MPa. J. Chem. Thermodyn.

2002, 34, (6), 717-734. 51.

Claus, P.; Kleinrahm, R.; Wagner, W., Measurements of the (p,ρ,T) relation of ethylene, ethane, and sulphur hexafluoride in the temperature range from 235K to 520K at pressures up to 30MPa using an accurate single-sinker densimeter. J. Chem. Thermodyn. 2003, 35, (1), 159-175.

52.

Fatoorehchi, H.; Rach, R.; Tavakoli, O.; Abolghasemi, H., An Efficient Numerical Scheme to Solve a Quintic Equation of State for Supercritical Fluids. Chem. Eng. Commun. 2014, 202, (3), 402-407.

53.

Lee, J. H.; Kang, J. W.; You, S.-S.; Lee, C. S., An algorithm for convergent solutions in vapor–liquid equilibrium calculations using an equation of state. Fluid Phase Equilib.

2005, 234, (1–2), 34-41. 54.

Frenkel, M.; Chirico, R. D.; Diky, V.; Muzny, C. D.; Kazakov, A. F.; Magee, J. W.; Abdulagatov, I. M.; Kroenlein, K.; Diaz-Tovar, C.; Kang, J. W.; Gani, R. NIST ThermoData Engine, NIST Standard Reference Database 103b-Pure Compounds, Binary


Page 31 of 47

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


Mixtures, andChemical Reactions, 7.3; National Institute of Standards and Technology: 2012. 55.

Diky, V.; Chirico, R. D.; Muzny, C. D.; Kazakov, A. F.; Kroenlein, K.; Magee, J. W.; Abdulagatov, I.; Frenkel, M., ThermoData Engine (TDE): Software Implementation of the Dynamic Data Evaluation Concept. 9. Extensible Thermodynamic Constraints for Pure Compounds and New Model Developments. J. Chem. Inf. Model 2013, 53, (12), 34183430.

56.

Koak, N.; de Loos, T. W.; Heidemann, R. A., Effect of the Power Series Dispersion Term on the Pressure−Volume Behavior of Statistical Associating Fluid Theory. Ind. Eng. Chem. Res. 1999, 38, (4), 1718-1722.

57.

Polishuk, I., Addressing the issue of numerical pitfalls characteristic for SAFT EOS models. Fluid Phase Equilib. 2011, 301, (1), 123-129.

58.

Polishuk, I.; Mulero, A., The numerical challenges of SAFT EoS models. Rev. Chem. Eng.

2011, 27, (5-6), 241. 59.

Jackson, G.; Gubbins, K. E., Mixtures of associating spherical and chain molecules. Pure Appl. Chem. 1989, 61, (6), 1021-1026.



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 32 of 47

APPENDIX

The original two parameters vdW EOS3 can be represented as follows: P=

RT a − 2 v−b v

(A.1)

where P is the pressure, T is the absolute temperature, v is the molar volume of the fluid, and R is the universal gas constant. The constants a and b are determined by applying critical constraints demanding the critical isotherm to pass through the critical point with zero values of the first and second derivatives of the pressure with respect to volume

a=

27 RTc 64 Pc

bc = b(Tc ) =

(A.2)

RTc 8Pc

(A.3)

where subscript c denotes the critical point. If the volume fraction is defined by φ = b/v, at the critical temperature, the substitution of eqs A.2 and A.3 into eq A.1 leads to a reduced form, Pbc 1 (3ϕ − 1)3 1 = + RTc 8 (1 − ϕ ) 8

(A.4)

where the left-side term indicates the reduced pressure and bc is the hard-core volume at critical temperature. A derivation of eq A.4 can be made through a different methodology, which begins with the assumption of the desirable mathematical form of the given EOS at critical state. If the critical volume fraction is denoted by β, for the original vdW EOS, the reduced critical isotherm satisfying the critical constraints is expected to have the form of Pbc F (ϕ / β − 1)3 = +F RTc 1− ϕ

(A.5)


Page 33 of 47

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


where the physical constraint of zero pressure at zero density is also met. This proposed form should satisfy the following two constraints: Constraint 1 from the ideal gas behavior,

Pbc b ∂P = lim c =1 ϕ→0 RT ϕ ϕ→0 RT ∂ϕ c c

lim Z = lim ϕ→0

(A.6)

where the second equality is based on l’Hospital’s rule and the relation of F = β/(3-β) is derived. Another constraint, Constraint 2, can be obtained through a simple mathematical operation; if the repulsive term, φ/(1 − φ), is subtracted from both sides of eq A.5, 3 Pbc ϕ F (ϕ / β − 1) ϕ − = +F− RTc 1 − ϕ 1− ϕ 1− ϕ

F (ϕ / β − 1) + F (1 − ϕ ) − ϕ = 1− ϕ 3

(A.7)

We note that the right-side term of eq A.7 corresponds to an attractive term of vdw EOS at the critical temperature. As the attractive term has no singularity at φ = 1, the numerator of the rightside of eq A.7 should contain 1 − φ for a cancellation with the denominator, 1 − φ. This reduces to the following relation, F (1 / β − 1) − 1 = 0 3

(A.8)

If we substitute F = β/(3 − β) into eq A.8, it is rearranged into

 1 3 3  β 3 − 2 + − 1 = 3 − β β β β 

(A.9)

This equation is reduced to 1-3β = 0. With β = 1/3, we obtain F = 1/8 and eq A.4 is again recovered.



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 34 of 47

Table 1. Optimized model parameters of four different EOS forms

Model parameters EOS form

bc

AAD(%)

m1

m2

m3

c1

c2

k

Eq 4

1.9068

6.3

-0.99

-

-

-

0.5436

53.5938

2.2

Eq 4

2.0094a

1.0094a

-0.96

-

-

-

0.4884

48.1669

4.2

Eq 19

-1.7162

0.9836

-

-3.8648

3.7992

-

0.4302

42.4133

1.7

Eq 28

-1.9695

4.7030

-

-4.7121

7.5912

-0.3273

0.3491

34.4213

1.0

a

β [dm3/mol]

These values correspond to m1*= 1.0094 and m2* = 1.


Page 35 of 47

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


Table 2. The comparison of AAD of both models for nine components

BWRS EOS

Eq 35

Ref

Substance ψc

e4

AAD(%)

k0

m1

m2

AAD(%)

CH4

0.054

0.835

0.68

-0.337

-2.039

4.731

0.97

31, 35a

C 2 H6

0.052

0.843

0.97

-0.315

-2.336

5.210

1.32

39b

C 3 H8

0.051

0.805

0.61

-0.326

-1.980

4.791

1.17

31, 37c

C 2 H4

0.050

1.027

1.62

-0.264

-2.697

5.830

2.21

[40, 41]d

Ar

0.054

0.868

2.29

-0.372

-2.206

4.724

2.39

42e

N2

0.053

0.924

0.89

-0.355

-2.259

4.861

1.03

43f

CO2

0.053

0.645

1.02

-0.198

-2.390

5.888

1.59

44-46g

SO2

0.057

0.472

1.36

-0.200

-2.100

5.610

1.46

47-49h

SF6

0.049

1.040

0.64

-0.267

-2.937

6.061

1.34

50, 51i

a

PVT data at 190.555 K were taken

b


c


d

In both refs, PVT data at 282.35 K were taken

e


f


g

PVT data at 304.135 K in ref 44 , data at 304.144 K in ref 45, and data at 304.187 K in ref

46 h

PVT data at 430.587 K

i

PVT data at 318.723 K in ref 50 and data at 318.7 K in ref 51



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 36 of 47

Table 3. The optimized model parameters

Non-critical parametersa

critical parametersb

Substance m1

m2

κ1

κ2

ε1

ε2

δ

Γ0,c

Γ1,c×10

bc [cm3/mol]

CH4

-1.918

4.268

0.450

0.125

1.010

0.148

0.499

3.4232

-3.424

35.35

C2H6

-2.183

5.086

0.851

0.163

1.329

0.349

0.502

3.6241

-2.863

50.70

C3H8

-1.884

4.671

0.858

0.214

1.292

0.335

0.503

3.6736

-3.574

69.81

C2H4

-2.640

5.691

0.819

0.185

1.378

0.492

0.503

3.5861

-2.138

45.18

Ar

-2.189

4.615

0.328

0.207

0.453

0.355

0.540

3.3781

-2.781

26.77

N2

-2.302

4.801

0.438

0.230

0.733

0.411

0.631

3.3876

-2.576

31.83

CO2

-2.538

5.954

0.726

0.405

1.087

0.589

0.610

3.7907

-2.362

31.61

SO2

-2.110

5.597

0.942

0.424

0.837

0.659

0.760

3.9614

-3.151

40.59

SF6

-2.973

6.134

0.832

0.291

1.221

0.474

0.621

3.5699

-1.733

68.43

a

the parameters not determined by solving both critical conditions and Zccal = Zcexp

b

the parameters determined by solving both critical conditions and Zccal = Zcexp


Page 37 of 47

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


FIGURE CAPTIONS

Figure 1. Critical isotherms of methane reproduced by eq 4 with two sets of parameter; the solid line with m1 = 1.9068, m2 = 6.3 and m3* = - 0.99 and the dashed line with m1 = 2.0094, m2 = 1.0094 and m3 = -0.96.

Figure 2. Critical isotherms of methane reproduced by eqs 19 and 4

Figure 3. Critical isotherms of methane reproduced by eqs 19 and 28

Figure 4. Critical isotherms of propane reproduced by eq 35 and BWRS EOS

Figure 5. Critical isotherms of SF6 reproduced by eq 35 and BWRS EOS

Figure 6. Density dependence of the terms constituting eq 36 for methane at Tc. The red dashed line indicates the contribution by the first term, the orange dashed line that by the second term, and the red solid line by eq 36

Figure 7. Effect of temperature dependent energy parameters on the behavior of eq 42 and its constituting contribution for methane at Tr = 0.3 (a) Γ0 = Γ0,c / Tr and Γ1 = Γ1,c / Tr0.5 (b) Γ0 = Γ0,c / Tr and Γ1 = Γ1,c / Tr2.5 The blue line represent the repulsive contribution



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

Figure 8. Reporduced reduced isotherms of carbon dioxide at three temperatures. Blue line indicates the isotherm at T = Tc /2, green line at T = Tc /5 and red line at T = Tc /10.

Figure 9. (a) Saturated vapor pressure curve for methane reproduced by the present model; , experimental data generated by REFPROP (b) P-ρ diagram for methane reproduced by the present model; , experimental data generated by REFPROP . Temperatures of the generated isotherms are 100, 120, 140, 160, 180, 190.56, 200, 220, 240, 260, 280 and 300K


Page 38 of 47

Page 39 of 47

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


Figure 1.



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

Figure 2.


Page 40 of 47

Page 41 of 47

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


Figure 3.



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

Figure 4.


Page 42 of 47

Page 43 of 47

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


Figure 5.



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

Figure 6.


Page 44 of 47

Page 45 of 47

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


Figure 7.



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

Figure 8.


Page 46 of 47

Page 47 of 47

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


Figure 9..


Development of an Accurate Quintic Equation of State: The Necessity

Recommend Documents