Nonphysical Behavior in Several Statistical Mechanically Based

2 days ago - Two studies published earlier this month, one from the biotech start-up Yumanity Therapeutics and another... BUSINESS CONCENTRATES ...
0 downloads 0 Views 2MB Size
Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Nonphysical Behavior in Several Statistical Mechanically Based Equations of State Nayef M. Alsaifi,*,† Mohammed Alkhater,† Housam Binous,‡ Isa Al Aslani,§ Yousef Alsunni,† and Zhen-Gang Wang⊥ †

Chemical Engineering Department, King Fahd University of Petroleum & Minerals, Dhahran 31261, Saudi Arabia Chemical Engineering Department, National Institute of Applied Sciences and Technology, University of Carthage, Carthage 1054, Tunisia § Saudi Aramco R&D Center, King Abdullah University of Techology (KAUST), Thuwal 23955-6900, Saudi Arabia ⊥ Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, United States

Ind. Eng. Chem. Res. Downloaded from pubs.acs.org by UNITED ARAB EMIRATES UNIV on 01/11/19. For personal use only.



ABSTRACT: In this work, we utilize bifurcation diagrams to study the role of mathematical artifacts in deteriorating the physical behavior in statistical mechanically based equations of state of pure fluids. We study the impact of common empirical approximations usually employed to overcome some of the mathematical and physical challenges such as the parametrization of mean field models or pair correlations functions at contact. The proposed diagrams elucidate how the reduced molar volume bifurcates with the variation of temperature at constant pressure. We generate bifurcation diagrams for the modified van der Waals equation of state (EOS) of Poole et al, SAFT-VR Mie, Soft-SAFT, CK-SAFT, and the original SAFT EOSs for spherical and nonspherical molecules. We find that the bifurcation diagram can serve as a useful schematic tool to reveal the unphysical PVT behavior, demonstrate the existence of physical and spurious two-phase separation regions, and illustrate how the number of molar volume roots vary with temperatures. Our method shows that the presence of unphysical branches can cause spurious two-phase separation regions and create erroneous behavior in the stability limit of vapor−liquid equilibrium. We demonstrate that the existence of customary and spurious phase envelopes is accompanied by S-shaped behavior in the volume-temperature bifurcation diagrams. The study reveals that none of the SAFT models is free from producing unphysical behavior. While the SAFT-VR Mie EOS exhibits solid−liquid-like behavior for nonspherical molecules, the CK-SAFT EOS shows liquid−liquid demixing behavior for spherical and nonspherical compounds. For the soft-SAFT EOS, three different two-phase separation regions are observed in addition to the common vapor−liquid phase separation region.

1. INTRODUCTION

functional forms inevitable. For instance, various empirical analytical expressions have been proposed for pair correlation functions at contact with adjustable parameters obtained with the aid of integral equation theories or molecular simulation.16−19 The difficulty in evaluating the relevant integrals has also produced approximate expressions lacking obvious physical intuition.20−22 In addition, significant effort has alternatively been devoted to the usage of available accurate molecular simulation data in constructing empirical expressions describing PVT and other thermodynamic properties of various molecular interactions.23−26 In other cases, some theoretical models, which have been created to describe specific molecular interactions, were combined with empirical or semiempirical equations of state.27−29

The development of simple thermodynamic models has been a subject of great interest since the remarkable work of van der Waals.1 Much of the early effort was primarily dedicated to proposing simple empirical2,3 or semiempirical models.4,5 Because of their limitations and lack of accuracy, the interest has dramatically shifted toward statistical mechanics from which more involved analytical equations of state6−8 have been obtained using hard-sphere theories,9 perturbation methods,10,11 and integral equation theories.12,13 Even though first-principle approaches have generally been considered, empiricism has continued to play a significant role in developing theoretical models. Much of the empiricism is unfortunately hard to avoid due to the nonexistence of exact solutions for many statistical mechanical theories. In other cases, even if it is possible to obtain the exact solution, the derived mathematical expression is too involved to be used in practical applications.14,15 These drawbacks have made the creation of arbitrary and empirical © XXXX American Chemical Society

Received: September 23, 2018 Revised: December 31, 2018 Accepted: January 4, 2019

A

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. (a) Correct physical behavior (S-shaped curve) shown on the bifurcation diagram of reduced volume−temperature at subcritical reduced pressure (Pr = 0.4) from any equation of state without defects. The branch has liquid and gas regions where the two turning points bounding a mechanically unstable region. (b) Physical branches at different isobars (Pr = 0.4, 0.6, 0.8, 1, and 1.1). At the critical isobar (Pr = 1), there is only one turning point while there is no turning point above the critical isobar (Pr = 1.1). (c) PV spinodal and binodal curves for a typical pure compound from any equation of state where the dots on the spinodal curve representing the turning points in panel b.

temperature varies at different isobars to investigate the variation of multiple volume roots and show the existence of spurious two-phase separation regions in the perturbed-chain statistical association fluid theory (PC-SAFT).39 In the present work, we apply bifurcation diagrams to evaluate the impact of empirical approximations on the physical behavior of several statistical mechanically based equations of state. Our proposed volume-temperature bifurcation diagram is utilized to discover spurious two-phase separation regions, show the deterioration of the physical PVT behavior, demonstrate all unphysical branches and their impact, and show the variation of multiple volume roots with temperature. We elucidate how the combination of the physical and unphysical branches can cause the disappearance of physical spinodal points. Furthermore, we show that the unphysical branches can cause spurious two-phase separation regions. To cover various types of commonly used empirical approximations, we consider five different equations of state. We start with the modified van der Waals of Poole et al.27 to examine the combination of a simple semiempirical cubic equation of state and a statistical mechanical model of hydrogen bonds. We then adopt Wertheim’s first order thermodynamic perturbation theory40−43 (TPT1), which has been applied extensively and extended from simple systems44−48 to complex polymers.49,50 In this study, we focus on the most popular extension of the so-called statistical association fluid theory (SAFT) to evaluate the impact of the employed empirical approximations. We study SAFT-variable range based on the generic Mie interactions (SAFT-VR Mie),17,51 Soft-SAFT,52,53 the SAFT of Huang and Radosz54 with the mean field model of Chen and Kreglewski55 (CK-SAFT) and the SAFT of Chapman et al.8 (original SAFT). We create bifurcation diagrams to explore the unphysical behavior associated with these equations for both spherical and nonspherical pure fluids. We also compare their bifurcation diagrams with that of the PC-SAFT EOS.34,39 For more details about the selected models, the reader may consult their respective references. We briefly summarize their mathematical formulations in Appendix A.

Despite the success gained from applying statistical mechanically based equations of state, users have reported abnormal behaviors, which were not previously observed in the widely used simple semiempirical models. Among these are prediction of unrealistic thermo-physical properties,30,31 multiple unphysical volume roots,32−36 artificial two-phase separation regions,30,34,37,38 and the intersection of isotherms.31 The unphysical behavior reported with some theory-based models exists either within or beyond the practical region of applications. While the appearance of any pathological behavior within the practical region is unquestionably considered as a failure of the model, it is surprising that the existence of unphysical behavior outside the practical region influences the model accuracy within the practical region. For instance, Yelash et al.37,38 elaborated with examples that the existence of the spurious liquid−liquid coexistence regions in the PC-SAFT model, which was located outside fluid region at low temperature and high density, deteriorated the PVT accuracy of the model compared to experiment within the fluid region. Although this conclusion may not be generalized to other models, the influence of mathematical artifacts remains unclear and not carefully assessed in statistical mechanically based equations of state. A more important issue resides on the methodology needed to identify the defects caused by inherent mathematical artifacts. In previous studies,30,31,33,36−38 the discovery of model defects was apparently found in an ad-hoc manner and no detailed methodology was proposed for this purpose. Indeed, it becomes a laborious task to find the deficiency in a given model by studying random values of state variables such as temperature, pressure, and density. Furthermore, the role of the proposed empirical approximations cannot be evaluated in a detailed and rigorous way. A more effective and appropriate approach is necessary to examine theoretical models and guarantee the safety of their applications under any state conditions. An attractive alternative lies in the construction of bifurcation diagrams, which enable us to study in detail the physical and unphysical behavior associated with any thermodynamic model. In a recent study,34 we proposed bifurcation diagrams in which molar volume bifurcates as B

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

corresponding two-phase coexistence. This can be verified from Figure 1b and c and the discussion in the next section. In Figure 1a, we only show one S-shaped curve, which represents one physical fluid−fluid transition (vapor−liquid phase transition). This should not imply that some physical equations of state cannot support more than one physical transition for pure fluids. In fact, the modified van der Waals of Poole et al.,27 which will be studied in Section 4, was constructed to support two physical fluid−fluid transitions for water. Nonetheless, the conclusions drawn from the physical behavior of water cannot be generalized to other pure fluids. There is no physical evidence to support the existence of more than one physical fluid−fluid transition in other pure compounds. Therefore, most equations of state including the SAFT EOSs are supposed to produce only the vapor−liquid phase transition for pure fluids. In addition, these equations were developed to describe fluid regions. We thus consider any additional phase transition to be spurious.

We organize the work as follows. Section 2 outlines the general features of the bifurcation diagram from any equation of state that is free from any defect. In Section 3, we give detailed analysis on how to detect multiple fluid-phase coexistence of pure compounds. Sections 4 and 5 explore the bifurcation diagrams of the selected equations of state. Discussion and Conclusions are given in the last two sections of the paper.

2. GENERAL FEATURES OF BIFURCATION DIAGRAM AND ITS APPLICATIONS We begin with introducing the general features of bifurcation diagrams generated from equation of state by pseudo arclength continuation method.34,56,57 The reader may consult our previous work where we described the details of the numerical analysis.34 A brief summary is given in Appendix B. Our proposed bifurcation diagram gives a detailed description of every possible volume root obtained from equation of state as the temperature varies at a fixed pressure. Figure 1a is a qualitative demonstration of the bifurcation diagram at a subcritical pressure (Pr = 0.4), where the axis of ordinates represents reduced molar volume and the axis of abscissas corresponds to reduced temperature. If the model is free from any physical defect, the only branch generated at any subcritical pressure is the S-shaped curve depicted in Figure 1a. We formerly confirmed this behavior by a comparison with experiment and analyzing the expected physical behavior of pure fluids over a wide range of state conditions.34 The S-shaped branch contains the loci of gas, liquid, and unstable volume roots at a wide range of temperatures at a fixed pressure. For instance, the bifurcation diagram (Figure 1a) shows only one liquid volume root at Tr = 0.6, while it gives the common three volume roots at Tr = 0.8. The two dots on the branch are called turning points and they bound the mechanically unstable region. The behavior of other isobars is shown on Figure 1b. The two turning points approach each other as the pressure increases, and they eventually collapse at the critical isobar (Pr = 1). Above the critical point, the branch has a sigmoid function shape with no turning point as seen for the branch at Pr = 1.1. To get more physical insight about the turning points, we show the binodal and spinodal curves in Figure 1c, where the points on the spinodal curve represent the turning points given in Figure 1b. To follow the turning points easily, we number them from 1 to 7 on both Figure 1b and c. The advantage of generating the bifurcation diagrams is not limited to the representation of volume roots. In fact, it provides unambiguous evidence of the existence of defects in any proposed model. This is because the unphysical behavior in any model arises either from the existence of additional branches other than the physical one or the extension of the physical branch beyond the S-shaped curve. For instance, if an unphysical branch exists close to the physical one, a change of pressure may lead to the combination of both branches. The combination may deteriorate the PVT behavior, form spurious two-phase separation regions, and delete physical spinodal points as will be demonstrated in Section 5. Another advantage of the bifurcation diagram is that it provides evidence of the existence of two-phase separation regions of pure compounds without conducting equilibrium calculations. As long as there is an S-shaped branch (or sometimes a flipped S-shaped curve) with two turning points separating the mechanically stable and unstable regions in the way shown in Figure 1a, there will be a

3. MULTIPLE FLUID-PHASE COEXISTENCE OF PURE COMPOUNDS The existence of spurious two-phase separation regions, besides the common vapor−liquid phase separation region, was reported in theory-based models.30,33,34,37,38 It is surprising that these spurious two-phase regions share the common thermodynamic behavior found in the physical fluid-phase coexistence region. Here, we give some guidelines, to be employed later, on how to pinpoint the artificial two-phase separation regions, if any, in any theory-based model. The guidelines are constructed based on our results in the present study and in our previous paper.34 Unlike the determination of the physical vapor−liquid coexistence region, it is not an easy task to confirm the existence of the spurious two-phase separation region in theory-based equations of state. The reasons are that such a behavior is not expected in a theorybased model and there is also no clue at what state conditions under which such phenomena may arise. This explains why it took years to clearly identify this problem in some theoretical equations of state.37 In the subsequent sections, we report this problem for the selected models; one of these models was previously thought to be free of such defect.17 An appropriate methodology of how to discover the existence of spurious two-phase regions has not been reported in the literature. There are, however, obvious common methods that could be utilized for this purpose with varying levels of difficulty. One such approach is to search for the critical points by solving algebraic system of equations satisfying the critical conditions. This approach is not applicable if the model has solid−liquid-like behavior where no critical point exists. An example of this phenomenon will be given later in Section 5.1. An alternative and simpler method is to test a random isotherm on the pressure−volume diagram and then inspect the van der Waals loop. Another option is to carry out equilibrium calculations at random state conditions. Although the latter two methods might work, randomly selecting of temperature is not efficacious. The reason is that the spurious two-phase envelope sometimes spans over a narrow temperature range, which can be easily missed by using arbitrary guesses. This is actually our motivation of using temperature as a bifurcation parameter. We thus propose to use the volumetemperature bifurcation diagram to observe any potential twophase transitions at a specific isobar. As previously indicated, C

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research the appearance of the S-shaped curve (or sometimes a flipped S-shaped curve) is sufficient to confirm its existence. The position of the S-shape curve of the spurious two-phase regions is most commonly found along with the extension of the physical branch, although it may be formed in a separate and unphysical branch. Therefore, it should not be surprising to observe two S-shaped curves on a single branch. Furthermore, we find from all studied models that the spurious two-phase separation regions expand over a broad range of pressure values. This explains why we propose bifurcation diagrams at a fixed isobar and not at a specific isotherm. Finally, the combination of unphysical and physical branches may cause the formation of spurious two-phase separation regions if certain conditions are satisfied. To avoid confusion, we postpone this subject to Section 5.2 where this formation is shown by an example.

Figure 2. Bifurcation diagram of reduced volume-temperature for water using the modified van der Waals EOS of Poole et al. at Pr = 0.05. The solid and dashed lines represent the mechanically stable and unstable volume roots; respectively. The flipped S-shaped curve magnified in the circle represents the liquid−liquid phase separation region. The adjustable parameters used for this figure are given by Poole et al.27 where εHB = −22 kJ/mol.

4. MODIFIED VAN DER WAALS EOS OF POOLE ET AL. Our first studied model is an example of the combination of a statistical−mechanical model of hydrogen bonds with a simple semiempirical equation of state. The model was proposed to study the anomalous properties of water by developing simple theory to describe the behavior of hydrogen bonds. The theory was based on partitioning the system into weak hydrogen bonds with zero energy and strong hydrogen bonds with εHB energy. It was then assumed that there are Ω ≫ 1 configurations of a weak bond and only one configuration of a strong bond. The Helmholtz free energy of hydrogen bonds was given by

also give the adjustable parameters of ethane to be used with the Soft-SAFT EOS in Section 5.2 to show the erroneous behavior caused by the combination of physical and physical branches. In the bifurcation diagrams, we scale the molar volume and temperature by critical constants obtained from experiment. 5.1. SAFT-VR Mie EOS. The SAFT-VR Mie EOS is one of the most recent versions of the SAFT family. Although there are various SAFT-VR models, this version is perhaps the most accurate and yet the most involved. The thermodynamic properties obtained from this model are in excellent agreement with molecular simulation data for monomer and chain fluids. The model showed improvement over other SAFT versions in predicting more accurate thermodynamic second-order derivative properties and phase equilibria.17,59,60 Nevertheless, the use of the Mie potential as a reference in the perturbation theory necessitated the use of approximate algebraic expressions for pair-correlation function at contact and for the mean field term. Although it was reported that the model is free from physical defects such as spurious two-phase separation region,17 there has been no rigorous demonstration to substantiate this claim. We generate reduced volume−temperature bifurcation diagrams for methane representing spherical molecules (no chain connectivity term) and n-octane representing nonspherical molecules at atmospheric pressure. For spherical particles, the upper branch in Figure 3 gives the correct physical behavior with two turning points without erroneous behavior. The five lower branches are unphysical, and some of them are mechanically stable (solid curves). Thus, the profusion of branches leads to multiplicity of the number of volume roots, whose number could reach a maximum of ten. Nonetheless, this is not problematic because their corresponding volume roots are located at low reduced molar volumes far from the fluid region. For other spherical molecules, it is not necessary to have the same number of lower branches. We find that the number of lower branches depends on the values of repulsive and attractive exponents (λr and λa) of the Mie potential.

−β AHB = f ln[Ω + exp( −βεHB)] + (1 − f) ln[Ω + 1] (1)

where f is the fraction of strong hydrogen bonds and it depends on density. This term was combined with the van der Waal equation of state, and the resulting equation is called the modified van der Waals equation of state of Poole et al.27 This equation displays a liquid−liquid phase envelope at low temperatures, in addition to the customary vapor−liquid equilibrium. Our volume−temperature bifurcation diagram illustrates the existence of this phase transition as well as the vapor−liquid phase transition on the same branch as shown in Figure 2. The liquid−liquid phase transition (the flipped Sshaped curve in the figure) located at about Tr = 0.2 is magnified and inserted in the same figure. While the used isobar to construct Figure 2 is at Pr = 0.05, the flipped Sshaped curve also appears at lower and higher pressures even above Pr = 1. We point out that the fraction of strong hydrogen bonds (f), which appears in eq 1, is an exponential expression of density given without rigorous definition. Nonetheless, we cannot confirm whether the liquid−liquid phase transition in water is physical or not since this issue is debatable.58 However, our bifurcation diagram clearly displays any two-phase separation region.

5. BIFURCATION DIAGRAMS OF VARIOUS SAFT MODELS We construct the reduced volume-temperature bifurcation diagrams of SAFT-VR Mie, Soft-SAFT, CK-SAFT, and original SAFT EOSs for spherical and nonspherical components such as methane and n-octane, respectively. The adjustable parameters of these compounds are given in Table 1. We D

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Table 1. Adjustable Parameters of Methane, Ethane, and n-Octane for SAFT-VR Mie, Soft-SAFT, CK-SAFT, and Original SAFT compound

model

m

σ (Å)

ε/k (K)

λr

λa

ref

methane

SAFT-VR Mie Soft-SAFT CK-SAFT Original SAFT Soft-SAFT SAFT-VR Mie Soft-SAFT CK-SAFT Original SAFT

1.0000 1.0000 1.0000 1.0000 1.5936 2.6253 3.9476 6.0450 8.0000

3.7412 3.7220 3.7004 3.7390 3.5850 4.4696 3.7870 3.0627 2.9766

153.360 147.302 190.290 152.200 147.302 369.180 249.811 206.030 148.100

12.650

6

17.378

6

17 52 54 8 52 17 52 54 8

ethane n-octane

Figure 4. Bifurcation diagram of reduced volume−temperature for noctane using the SAFT-VR Mie EOS at atmospheric pressure. The solid and dashed lines represent the mechanically stable and unstable volume roots, respectively.

reported in other models like the PC-SAFT EOS, there is a major distinction that makes it unique. Unlike other models displaying liquid−liquid-like transition, the SAFT-VR Mie displays a solid−liquid-like transition and no critical point is found although we test the model at very high pressure. The spurious phase envelope of n-octane is shown in Figure 5b in addition to the common vapor−liquid envelope. In the same figure, we see four points (1−4) on the spinodal curve. Their corresponding turning points are demonstrated in the flipped S-shaped curves in Figure 5a at Pr = 0.1 and 140. The S-shaped curves, which belong to the spurious two-phase region, can be observed at any pressure in the V−T bifurcation diagram. On the other hand, if one adopts equilibrium calculations, the only possible region that could be determined is Region I depicted in Figure 5b. Outside this region, the values of fugacity are either extremely very big or almost zero. The spurious saturation pressure corresponding to Region I is illustrated in Figure 6 where the axis of ordinates corresponds to reduced saturation pressure, the axis of abscissas represents reduced temperature and the dashed lines give the spinodal curves. We believe the spurious two-phase region found in SAFTVR Mie is a result of a mathematical artifact of the pair correlation function at contact employed in the chain connectivity term since this phenomenon is only observed for chain-like molecules.

Figure 3. Bifurcation diagram of reduced volume−temperature for methane using the SAFT-VR Mie EOS at atmospheric pressure. The solid and dashed lines represent the mechanically stable and unstable volume roots, respectively.

For nonspherical molecules, the chain connectivity deteriorates the behavior of the physical branch. A bifurcation diagram of reduced volume against reduced temperature is depicted in Figure 4 for n-octane. The number of the unphysical branches for n-octane is less than that of methane but share small reduced volumes exceeding the limit of the fluid region. The number of volume roots reaches a maximum of five. At low temperatures, the physical branch has two additional turning points bounding a mechanically unstable region. Its Sshaped curve confirms the existence of additional two-phase separation region in addition to the physical one. In consequence, the SAFT-VR Mie displays spurious two-phase separation region in contrast to what was previously advocated. The existence of the spurious two-phase separation region is not restricted to n-octane, but rather a general feature arising from the application of the SAFT-VR Mie to any nonspherical molecules. Although this phenomenon resembles what was E

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 5. (a) Two S-shaped curves in the V−T bifurcation diagram generated from SAFT-VR Mie theory for n-octane. (b) Physical and spurious phase envelopes for n-octane. The dots numbered from 1 to 4 represent the turning points in panel a. Region I is the selected region of the spurious saturation pressure depicted in Figure 6.

we investigate the impact of the employed multiparameter equations to soft-SAFT on physical behavior for both spherical and nonspherical molecules. For spherical molecules, the V−T bifurcation diagram of methane is given in Figure 7 at atmospheric pressure. The

Figure 6. Spurious phase diagram of n-octane generated from the SAFT-VR Mie EOS. The solid and dashed curves are saturation pressure and spinodal curves, respectively.

5.2. Soft-SAFT EOS. The Soft-SAFT EOS is a LennardJones (LJ) based model and it generally gives quite reasonable correlations for phase equilibria and thermophysical properties.61−63 A noticeable feature of the Soft-SAFT is the utilization of existing analytical models developed by accurate LJ simulation data. Among these models is the accurate approximation of the LJ model proposed by Johnson et al.,64 which constitutes the excluded volume and mean field term in the Soft-SAFT. The EOS of Johnson et al., which is a modified version of the Benedict-Webb-Rubin (MBWR) EOS,2 showed remarkable accuracy in correlating PVT and internal energies from the triple point up to 4.5-times the critical temperature. In fact, for spherical particles, the Soft-SAFT is mathematically identical to that of Johnson et al. The mathematical formulation of the EOS of the Johnson et al. is an involved functional form of density and temperature with 32 linear and one nonlinear parameters. To account for monomers bonding, an additional empirical equation was utilized for the LJ cavity correlation function at contact with 25 adjustable parameters.16 Some of the defects associated with multiparameter models were reported in the literature for various models.65,66 Here,

Figure 7. Bifurcation diagram of reduced volume−temperature for methane using the Soft-SAFT EOS at atmospheric pressure. The solid and dashed lines represent the mechanically stable and unstable volume roots, respectively.

physical branch shares the same defect observed also in the SAFT-VR Mie for nonspherical particles where two additional turning points exist at low temperatures. Hence, the existence of spurious two-phase separation region is evident as a result of the flipped S-shaped curve. The Soft-SAFT also accommodates unphysical isola (a closed and isolated branch) with two turning points located below the physical branch as seen in the figure. While the upper portion of the isola is mechanically unstable, the lower part is stable. At about reduced temperature of 0.5, we place two adjacent points on both branches to emphasize the importance of exercising caution in the calculations of volume roots. It is interesting to note that this problem is observed above the triple point of methane (Tr ≈ 0.47) where the EOS of Johnson et al. is quite accurate in F

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

turning point denoted by bifurcation point in Figure 9b. A further decrease of the pressure beyond this point shows erroneous behavior of the physical branch and no spinodal limit of the vapor phase is observed as shown in Figure 9c at Pr = 0.02. It is important to note that the corresponding pressure of Pr = 0.02 is equivalent to 1 atm. For the purpose of clarity of the expose, we do not show the lower branches of the V−T bifurcation diagrams in Figure 9. We also illustrate another problem for ethane in Figure 10 where the combination of unphysical and physical branches leads to the appearance of spurious two-phase separation regions. In Figure 10a, the V−T bifurcation diagram is similar to that of Figure 7 where we notice an isola below the physical branch. The increase of pressure forms an independent upper branch with four turning points while the physical branch moves down toward the isola (Figure 10b). The separation between the two branches is shown in the zoomed-in part. At Pr = 75.9, the physical branch touches the unstable portion of the isola and forms a bifurcation point as shown in Figure 10c. Any further increase of the pressure leads to the formation of an S-shaped curve on the right and a flipped S-shaped curve on left (Figure 10d). These two S-shaped curves prove the existence of new spurious two-phase separation regions. In Figure 10a and d, we observe four S-shaped curves numbered from 1 to 4. The first S-shaped curve is physical, while the other three are unphysical. Therefore, the Soft-SAFT EOS displays up to three spurious coexistence phase envelops. They are not shown here but we find that one of the spurious two-phase separation region is solid−liquid-like and the other two are liquid−liquid-like. The spinodal points on the third flipped S-shaped curve, shown on Figure 10d, belong to the solid−liquid-like phase separation region. The formation of spurious two-phase regions, as a result of the combination of the physical and unphysical branches, occurs not only for ethane or n-octane but also for methane. In Figure 7, we have found one spurious two-phase separation region due to the formation of the S-shaped curve at low temperature and 1 atm. However, if the pressure is increased, the physical branch touches the unphysical one in the same way demonstrated in Figure 10. As a result, two additional Sshaped curves are observed for methane. Therefore, the problem of the existence of the unphysical branches should not be underestimated, particularly if these branches are close to

reproducing simulation data. It is also important to realize that the mechanical stability test is not sufficient due to the existence of stable unphysical branches. For nonspherical molecules, the addition of the chain connectivity causes deterioration of the physical PVT behavior at low temperatures as shown in Figure 8 for n-octane at

Figure 8. Bifurcation diagram of reduced volume−temperature for noctane using the Soft-SAFT EOS at atmospheric pressure. The solid and dashed lines represent the mechanically stable and unstable volume roots, respectively.

atmospheric pressure. The bifurcation diagram shows 12 turning points and a maximum of 8 volume roots. The complexity of unphysical branches is illustrated in the zoomedin part in the figure. The unphysical behavior, due to the inclusion of the chain connectivity term, is not limited to noctane but also seen in short-chain molecules such as ethane. In fact, there are more serious problems in shorter chain molecules due to the combination of unphysical branches with the physical one. We demonstrate one of the problems using the V−T bifurcation diagrams of ethane in Figure 9. The unphysical branch on the left (Figure 9a) at Pr = 0.1 does not remain in its position if the pressure is decreased. At Pr = 0.03, this unphysical branch touches the physical one at the upper

Figure 9. Bifurcation diagrams of reduced volume−temperature for ethane using the Soft-SAFT EOS at (a) Pr = 0.1, (b) Pr = 0.03, and (c) Pr = 0.02. The solid and dashed lines, respectively, represent the mechanically stable and unstable volume roots. The figures show how the combination of physical and unphysical branches deteriorates the physical PVT behavior when pressure changes. The combination of the two branches deletes the vapor spinodal point at 1 atm. G

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 10. Bifurcation diagrams of reduced volume-temperature for ethane using the Soft-SAFT EOS at (a) Pr = 0.04, (b) Pr = 13, (c) Pr = 75.9, and (d) Pr = 102. The solid and dashed lines respectively represent the mechanically stable and unstable volume roots. The figures show how the combination of physical and unphysical branches causes the formation of spurious two-phase separation regions.

Figure 11. Bifurcation diagram of reduced volume-temperature for n-octane using (a) CK-SAFT, (b) original SAFT, and (c) PC-SAFT at atmospheric pressure. The solid and dashed lines represent, respectively, the mechanically stable and unstable volume roots.

the physical one. This problem is not only limited to SoftSAFT but we also observed it in PC-SAFT.34 5.3. Original SAFT and CK-SAFT EOSs. The original SAFT and CK SAFT equations of state played a major role in the foundation of other SAFT models and they have extensively been applied to various types of mixtures.8,67,69 The major difference between these two models lies in the mean field term. In the CK-SAFT, the mean field term was empirically expanded as a power series in terms of density and temperature. The coefficients were adjusted to PVT data of argon. On the other hand, the mean field term in the original SAFT had a much simpler expansion of density and temperature with coefficients obtained from molecular simulation data of Lennard-Jones fluid. These two mean field

terms in the two models represent their major source of empiricism. We generate V−T bifurcation diagrams from both models for n-octane at 1 atm (Figure 11a,b). For the sake of comparison, we add the bifurcation diagram of PC-SAFT (Figure 11c), which utilizes the same excluded volume and chain terms; however, it mathematically differs from the other two models in the mean field term. From the bifurcation diagrams, we see that the CK-SAFT and PC-SAFT display spurious two-phase separation regions due to the existence of the S-shaped behavior located on the physical branch at low temperatures. The original SAFT is surprisingly free from this defect. The spurious two-phase separation regions associated with these models were previously reported.30,31,34,37 However, H

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

demonstrate the superiority of the bifurcation analysis to commonly used previous methods. In Table 2, we give a summary of a list of physical defects by numbers in SAFT models as a result of the mathematical

our approach detects this defect without the need for any random testing of state variables. From the figures, we also observe the multiplicity of volume roots caused by the presence of unphysical branches. For spherical molecules, the CK-SAFT shares similar bifurcation diagram to Figure 11a. On the other hand, the PC-SAFT and original SAFT have a fewer number of unphysical branches for spherical particles. In the absence of the chain term, the branches designated by A and B in Figure 11b and c are absent.

Table 2. Summary of Unphysical Behavior in SAFT-VR Mie, Soft-SAFT, CK-SAFT, Original SAFT, and PC-SAFT EOSs for Spherical and Nonspherical Pure Fluids at 1 atm spherical pure fluids (without chain connectivity term) SAFT-VR Mie

6. DISCUSSION We utilize V−T bifurcation diagrams to detect the unphysical behavior in theoretical equations of state. We cover five equations of state with different types of empirical approximations. The modified van der Waals EOS of Poole et al. consists of the celebrated van der Waals equation of state with the addition of a statistical-mechanical model of hydrogen bonds. The other studied equations of state are SAFT-VR Mie, Soft-SAFT, PC-SAFT, CK-SAFT, and the original SAFT, which are all based on Wertheim’s first order thermodynamic perturbation theory. In each of these models, different empirical approximations have been used to overcome some of the mathematical or physical challenges. For the modified van der Waals EOS of Poole et al., our objective was mainly to demonstrate the S-shaped behavior associated with multiple phase separation regions in the V−T bifurcation diagrams since the model is quite simple. From a mathematical point of view, the original SAFT, CKSAFT, and PC-SAFT EOSs share the same excluded volume and chain connectivity models but differ substantially in the mean field terms regardless of which one is physically more realistic. In these three models, a large number of adjustable parameters are used in the mean field term: 42 in PC-SAFT, 24 in CK-SAFT, and 8 in the original SAFT. Irrespective of the accuracy, the original SAFT shows the correct PVT behavior without any spurious liquid−liquid separation region and it gives the correct number of turning points, although it still exhibits either one or two unphysical branches for spherical and nonspherical fluids, respectively. This indicates that the mathematical simplicity of the empirical approximation reduces the defects in PVT behavior. Besides model simplicity, the type of functional form is another important factor that influences the physical behavior of models. Although 47 coefficients were utilized in the correlation of the mean field term in SAFT-VR Mie, we find that the SAFT-VR Mie EOS gives the correct physical branch for spherical particles, albeit there are several unphysical branches at low reduced molar volumes. Simpler functional forms have been utilized in the CK-SAFT and PC-SAFT EOSs for the mean field term. Nevertheless, they contain spurious two-phase separation region for spherical fluids. The proposed functional forms were a rational expression in SAFT-VR Mie, polynomial expressions in CK-SAFT and PC-SAFT, and a combination of polynomial and exponential expressions in density in Soft-SAFT. For nonspherical particles where the chain connectivity term requires pair correlation function at contact, the bifurcation diagrams show more defects in the physical behavior in all studied models including SAFT-VR Mie, which displays spurious solid−liquid-like behavior of nonspherical molecules within a narrow temperature range. Our proof of the existence of this behavior in SAFT-VR Mie is in contrast to what has been concluded before.17 Our results

unphysical branches turning points spurious phase envelope unphysical critical points max volume roots

soft-SAFT PC-SAFT34

CK-SAFT

original SAFT

5

1

2

1

1

4

6

6

7

2

0

1 (3)a

1 (2)a

1

0

0

2

2

1

0

10

7

6

5

4

nonspherical pure fluids (with chain connectivity term) SAFT-VR Mie unphysical branches turning points spurious phase envelope unphysical critical points max volume roots

soft-SAFT PC-SAFT

CK-SAFT

original SAFT

2

5

3

1

2

4 1

12 1 (3)a

7 1

7 1

2 0

0

2

1

1

0

5

8

6

5

5

a

Number between brackets represents the maximum number of spurious phase envelopes found in these models by studying the pressure effect on the bifurcation diagram. The other numbers in this table are obtained at 1 atm. In this work, the pressure effect is only studied for Soft-SAFT.

artifacts for both spherical and nonspherical pure fluids at atmospheric pressure. We point out that these numbers may vary with pressure as shown in Section 5.2 for Soft-SAFT. Although it is practically possible to overcome these problems, it remains an obstacle for robust and fast computation since more time is needed to obtain and test unphysical volume roots caused by unphysical branches. The defects reported in Table 2 do not imply that a fundamental problem exist in the SAFT theory. As previously indicated, these defects are caused by the employed empirical approximations. Simpler mean field and segment terms may, but not necessarily, avoid these problems. For instance, the Elliott−Suresh−Donohue (ESD) EOS,70 which can be recast into the SAFT framework,71 is free from these defects. Figure 1272 shows the bifurcation diagram of ethane generated by ESD model at 1 atm. As seen, the ESD model gives the correct S-shaped curve without any additional branch. In a future study, we will extend our bifurcation analysis to evaluate other models. We will also perform bifurcation analysis on mixtures and on the effect of polar and associative terms. To sum up, the remedy of the physical defects associated with SAFT models is needed for several reasons. First, the SAFT model is a perturbation theory which continuously I

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

We found that the unphysical branches create unphysical twophase separation regions, cause the disappearance of physical spinodal points, deteriorate the physical PVT behavior, and cause multiplicity in the number of volume roots. While the mechanical stability testing is crucial in determining volume roots, it is not sufficient because some of the unphysical branches are mechanically stable.



APPENDIX A. SUMMARY OF GOVERNING EQUATIONS OF SAFT MODELS The extended TPT1 can be given as follows: a RES = a MONO + a CHAIN

(A.1)

where a is reduced Helmholtz free energy, RES corresponds to the residual Helmholtz free energy, CHAIN corresponds to chain connectivity between segments, and MONO represents monomer segment, which is the contribution of the excluded volume (HS) and mean field term. For SAFT-VR Mie, PCSAFT, CK-SAFT, and original SAFT, the excluded volume was given by Carnahan−Starling:9

Figure 12. Bifurcation diagram of reduced volume-temperature for ethane using the ESD equation of state at atmospheric pressure. The solid and dashed lines represent respectively the mechanically stable and unstable volume roots.

accommodates new physical effects (e.g., electrolyte and polar) as perturbation terms, which may cause artifacts since the employed models most likely add more empirical approximations. Second, the SAFT models were supposed to be more physical and safer to extrapolate. However, our work shows that they can possess many unphysical branches. Furthermore, the influence of these defects on the accuracy remains uncertain. Hence, it is necessary to resolve these issues to ensure the safety and robustness of statistical mechanically based equations of state. Finally, we emphasize the importance of bifurcation diagram as an efficient tool to study the available shortcomings in any thermodynamic model.

a HS = m

4η − 3η2 (1 − η)2

(A.2)

where η is packing fraction. In the extended TPT1, the chain connectivity term is given in terms of pair correlation function at contact for the specified reference: a chain = (1 − m)ln g(σ )

(A.3)

The PC-SAFT, CK-SAFT, and original SAFT share the same chain connectivity term where the pair correlation function at contact was derived from Carnahan−Starling: 1

g(σ ) =

7. CONCLUSIONS In this work, we studied the effect of mathematical artifacts on the physical behavior of various statistical−mechanical-based models including the SAFT-VR Mie, Soft-SAFT, CK-SAFT, the original SAFT, and the modified van der Waals of Poole et al. equations of state. To do so, we generated bifurcation diagrams where the molar volume bifurcates with the change of temperature at a certain pressure for either spherical or nonspherical pure fluids. We found that the mathematical artifacts led to the appearance of unphysical branches giving rise to multiplicity of volume roots in all the SAFT models. In addition, they are responsible for the extension of the main branch beyond the S-shaped physical behavior. The extension of the branch resulted in the existence of spurious liquid− liquid-like separation regions in CK-SAFT, Soft-SAFT, and PC-SAFT for spherical and nonspherical particles. For the SAFT-VR Mie, our study showed that the SAFT-VR Mie contains solid−liquid-like separation region of nonspherical pure fluids but not for spherical ones. The original SAFT, on the other hand, exhibits a maximum of two unphysical branches, but it does not contain any spurious phase separation regions. This was a result of the simplicity of the employed mean field term. We also found that the SoftSAFT displays a liquid−liquid separation region for both spherical and nonspherical molecules at atmospheric pressure. The Soft-SAFT accommodates two additional two-phase separation regions if the pressure effect is studied. Finally, the presence of unphysical branches in the statistical mechanically based equations of state should not be ignored.

1 − 2η (1 − η)3

(A.4)

The pair correlation function at contact for the SAFT-VR Mie is approximated by ij g (σ ) g (σ ) yz + (βε)2 2 zzzz g Mie(σ ) ≈ g 0(σ ) expjjjjβε 1 j g (σ ) g 0(σ ) z 0 k {

(A.5)

where g0(σ), g1(σ), and g2(σ) are lengthy expressions in terms of density and diameter, attractive, and repulsive exponents. For the Soft-SAFT, the pair correlation function at contact is LJ-based and given by23 5

g LJ(σ ) = 1 +

5

∑ ∑ aij(ρ*)i (T*) j − 1 (A.6)

i=1 j=1

where ρ* and T* reduced density and temperature, respectively. The coefficients (aij) can be found elsewhere.23 The rest of terms are defined independently since they are not similar. SAFT-VR Mie

iji σ y λr i σ y λa yz u Mie(r) = Cεjjjjjjj zzz − jjj zzz zzzz kr{ { kk r {

The SAFT-VR Mie is based on the Mie potential:

C= J

λr ijj λr yzz j z λr − λa jjk λa zz{

(A.7)

λa / λ r − λa

(A.8) DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research where λr and λa are repulsive and attractive exponents. The mean field in the SAFT-VR is given by a mf = m(β a1 + (β)2 a 2 + (β)3 a3)

where ai and bi are temperature-dependent functions, and Gi indicates exponential functions in terms of density. Although the above equation looks simple, the expansion of coefficients makes it quite involved and lengthy. The function has 32 adjustable parameter correlated to LJ molecular simulation data.

(A.9)

where k is Boltzmann constant, and T is temperature. a1, a2, and a3 are first-, second-, and third-order perturbation terms of Barker and Henderson (BH) perturbation theory.68 The firstand second-order perturbation terms are, respectively, given by

CK-SAFT

The Helmholtz free energy of the mean field in the CK-SAFT is given by55 ÄÅ ÉÑi ÄÅ ÉÑ j Å ε ÑÅηÑ a mf = m ∑ ∑ Dij ÅÅÅÅ ÑÑÑÑ ÅÅÅÅ ÑÑÑÑ ÅÇ kT ÑÖ ÅÇ τ ÑÖ

a1 = C[x 0 λa(a1S(η ; λa)) + B(η ; λa) − x 0 λr(a1S(η ; λr)) + B(η ; λr)]

(A.10)

i

and

where ε is well-depth, τ = 0.74048, and Dij are universal constants obtained by fitting PVT, second virial coefficients, and internal energy data of argon.

1 a 2 = KHS(1 + χ )εC2[x 0 2λa(a1S(η ; 2λa) + B(η ; 2λa)) 2 − 2x 0 λa + λr(a1S(η ; λa + λr) + B(η ; λa + λr)) + x 0 2λr(a1S(η ; 2λr) + B(η ; 2λr))]

Original SAFT

For the original SAFT, the Helmholtz free energy is an LJbased model fitted to molecular simulation data:69

(A.11)

In A.10 and A.11, χ is a correction factor, and (1 − η)4 1 + 4η + 4η2 − 4η3 + η 4

K HS =

a mf = m (A.12)

+ 10.285ρR 3 )

+ 15.904ρR 3 )

and

TR =

zyz z 3 zz{



(A.15)

(A.16)

f1(V(s), T(s)) = 0

i=1

aiρ*i + i

ij dV yz i dT y jj zz + jjj zzz = 1 ds k { k ds { 2

i=1

2

(B.3)

It should be noted that B.3 is a differential equation, while B.2 is an algebraic equation. Hence, one need an initial condition to solve this DAEs. To do so, one can write T(s = 0) = T0 where T0 is chosen arbitrarily. Then the root of the equation defined by f1(V(s = 0), T0) = 0 is computed using

6

∑ biGi

(B.2)

and

(A.17)

The soft-SAFT equation of state is a LJ-based model given by 8

(B.1)

where V is molar volume, T is temperature, and f1 is specific for each particular equation of state. The pseudo arc-length numerical technique introduces a new variable, s, called the arc-length segment then solving the following system of differential-algebraic equations:

Soft-SAFT



i 6 yz zzη ρR = jjj kπ 2 {

f1(V, T) = 0

where ϕi,n are 47 adjustable coefficients obtained using VLE data.17

ares =

kT ε

APPENDIX B. BRIEF SUMMARY OF PSEUDO ARC-LENGTH CONTINUATION METHOD USED TO GENERATE BIFURCATION DIAGRAMS For a fixed pressure, any equation of state is given by

and n=3 ∑n = 0 ϕi , nα n) ( fi (α) = i = 1, ..., 6 (∑nn==64 ϕi ,nα n−3)

(A.22)

and

(A.14)

where Iλ(λ) and Jλ(λ) are functions in terms of Mie exponents and segment diameter, ηeff(η;λ) is a fourth order polynomial in packing fraction, and the polynomial coefficients are expressions in terms of λ. The third-order perturbation term is an exponential function:

1 ji 1 − α = Cjjj j λa − 3 λr − k

(A.21)

2 a disp 02 = ρR ( − 1.9075 + 9.9724ρR − 22.216ρR

(A.13)

where

(A.20)

2 a disp 01 = ρR ( − 8.5959 − 4.5424ρR − 2.1268ρR

9η(1 + η) ji 1 − η /2 zy B(η ; λ) = 12ηεjjj I (λ ) − J (λ)zzz 3 λ j (1 − η)3 λ z 2(1 − η) k {

a3 = −ε 3f 4(α)η x 0 3 exp(f 5(α)η x 0 3 + f 6(α)η2 x 0 6)

disp ε R jij disp a 02 zyz jja 01 + zz k jk TR z{

where

where d is temperature-dependent segment diameter. The B and a1S are given in terms of packing fraction and Mie exponents:

i 1 yz 1 − ηeff (η ; λ)/2 zz a1S(η ; λ) = −12εηjjj 3 k λ − 3 { (1 − ηeff (η ; λ))

(A.19)

j

(A.18) K

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research



any nonlinear equation solving method. Hence, we find the second initial condition V(s = 0) = V0 and we can proceed with the solution of the DAEs using the built-in function of Mathematica, NDSolve, or any other solver. Finally, we generate the bifurcation diagram of V(s) versus T(s) at the fixed pressure.



Article

REFERENCES

(1) Van Der Waals, J. D.; Rowlinson, J. S. On the Continuity of the Gaseous and Liquid States; Courier Corporation, 2004. (2) Benedict, M.; Webb, G. B.; Rubin, L. C. An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures I. Methane, ethane, propane and n-butane. J. Chem. Phys. 1940, 8, 334−345. (3) Benedict, M.; Webb, G. B.; Rubin, L. C. An Empirical Equation for Thermodynamic Properties of Light Hydrocarbons and Their Mixtures II. Mixtures of Methane, Ethane, Propane, and n-Butane. J. Chem. Phys. 1942, 10, 747−758. (4) Redlich, O.; Kwong, J. N. On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions. Chem. Rev. 1949, 44, 233−244. (5) Peng, D.-Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (6) Beret, S.; Prausnitz, J. M. Perturbed hard-chain theory: An equation of state for fluids containing small or large molecules. AIChE J. 1975, 21, 1123−1132. (7) Dickman, R.; Hall, C. K. Equation of state for chain molecules: Continuous-space analog of Flory theory. J. Chem. Phys. 1986, 85, 4108−4115. (8) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New reference equation of state for associating liquids. Ind. Eng. Chem. Res. 1990, 29, 1709−1721. (9) Carnahan, N. F.; Starling, K. E. Equation of state for nonattracting rigid spheres. J. Chem. Phys. 1969, 51, 635−636. (10) Barker, J. A.; Henderson, D. Perturbation theory and equation of state for fluids: the square-well potential. J. Chem. Phys. 1967, 47, 2856−2861. (11) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237−5247. (12) Percus, J. K.; Yevick, G. J. Analysis of classical statistical mechanics by means of collective coordinates. Phys. Rev. 1958, 110, 1. (13) Lebowitz, J. L.; Percus, J. K. Mean spherical model for lattice gases with extended hard cores and continuum fluids. Phys. Rev. 1966, 144, 251. (14) Tang, Y.; Lu, B. C.-Y. Direct calculation of radial distribution function for hard-sphere chains. J. Chem. Phys. 1996, 105, 8262− 8265. (15) Chang, J.; Kim, H. Analytical expression for the correlation function of a hard sphere chain fluid. Mol. Phys. 1999, 96, 1789−1794. (16) Johnson, J. K.; Mueller, E. A.; Gubbins, K. E. Equation of state for Lennard-Jones chains. J. Phys. Chem. 1994, 98, 6413−6419. (17) Lafitte, T.; et al. Accurate statistical associating fluid theory for chain molecules formed from Mie segments. J. Chem. Phys. 2013, 139, 154504. (18) Goldman, S. An explicit equation for the radial distribution function of a dense Lennard-Jones fluid. J. Phys. Chem. 1979, 83, 3033−3037. (19) Matteoli, E.; Mansoori, G. A. A simple expression for radial distribution functions of pure fluids and mixtures. J. Chem. Phys. 1995, 103, 4672−4677. (20) Gubbins, K. E.; Twu, C. H. Thermodynamics of polyatomic fluid mixturesI theory. Chem. Eng. Sci. 1978, 33, 863−878. (21) Nicolas, J. J.; Gubbins, K. E.; Streett, W. B.; Tildesley, D. Equation of state for the Lennard-Jones fluid. Mol. Phys. 1979, 37, 1429−1454. (22) Luckas, M.; Lucas, K.; Deiters, U.; Gubbins, K. E. Integrals over pair-and triplet-correlation functions for the Lennard-Jones (12−6)fluid. Mol. Phys. 1986, 57, 241−253. (23) Johnson, J. K.; Mueller, E. A.; Gubbins, K. E. Equation of state for Lennard-Jones chains. J. Phys. Chem. 1994, 98, 6413−6419. (24) Kolafa, J.; Nezbeda, I. The Lennard-Jones fluid: An accurate analytic and theoretically-based equation of state. Fluid Phase Equilib. 1994, 100, 1−34. (25) Müller, A.; Winkelmann, J.; Fischer, J. Simulation studies on mixtures of dipolar with nonpolar linear molecules II. A mixing rule

AUTHOR INFORMATION

Corresponding Author

*Tel.: +966-138602194. Fax: + 966-3-860-4234. E-mail: alsaifi@kfupm.edu.sa. ORCID

Nayef M. Alsaifi: 0000-0003-3232-6411 Zhen-Gang Wang: 0000-0002-3361-6114 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors gratefully acknowledge the support of King Fahd University of Petroleum & Minerals through the grant to KFUPM-Caltech Research in Catalysis.



NOMENCLATURE reduced Helmholtz free energy first-order perturbation term of Barker and Henderson theory a2 second-order perturbation term of Barker and Henderson theory a3 third-order perturbation term of Barker and Henderson theory ai Functions of temperatures in Soft-SAFT aij Coefficients in the LJ pair correlation function at contact used in Soft-SAFT bi Functions of temperatures in Soft-SAFT d temperature-dependent segment diameter Dij Universal constants in CK-SAFT g(σ) pair correlation function at contact Gi Functions of density in Soft-SAFT k Boltzmann constant m number of segments mf mean field T temperature T* reduced temperature a a1



GREEK SYMBOLS η packing fraction ϕi,n adjustable parameters in the third-order perturbation theory of BH ε well-depth σ diameter of segment λr repulsive exponent in Mie potential λa attractive exponent in Mie potential χ correction factor in the mean field term in SAFT-VR Mie ρ* Reduced density



SUB, SUPERSCRIPTS, AND ABBREVIATIONS BH Barker and Henderson CHAIN chain connectivity HS hard sphere MONO monomer segment RES residual Helmholtz free energy SAFT Statistical association fluid theory L

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research for the dipolar contribution to the Helmholtz energy. Fluid Phase Equilib. 1996, 120, 107−119. (26) Saager, B.; Fischer, J. Construction and application of physically based equations of state: Part II. The dipolar and quadrupolar contributions to the Helmholtz energy. Fluid Phase Equilib. 1992, 72, 67−88. (27) Poole, P. H.; Sciortino, F.; Grande, T.; Stanley, H. E.; Angell, C. A. Effect of hydrogen bonds on the thermodynamic behavior of liquid water. Phys. Rev. Lett. 1994, 73, 1632. (28) Saager, B.; Hennenberg, R.; Fischer, J. Construction and application of physically based equations of state: Part I. Modification of the BACK equation. Fluid Phase Equilib. 1992, 72, 41−66. (29) Kontogeorgis, G. M.; Voutsas, E. C.; Yakoumis, I. V.; Tassios, D. P. An equation of state for associating fluids. Ind. Eng. Chem. Res. 1996, 35, 4310−4318. (30) Polishuk, I. About the numerical pitfalls characteristic for SAFT EOS models. Fluid Phase Equilib. 2010, 298, 67−74. (31) Polishuk, I. Addressing the issue of numerical pitfalls characteristic for SAFT EOS models. Fluid Phase Equilib. 2011, 301, 123−129. (32) Alsaifi, N. M.; Englezos, P. Prediction of multiphase equilibrium using the PC-SAFT equation of state and simultaneous testing of phase stability. Fluid Phase Equilib. 2011, 302, 169−178. (33) Privat, R.; Gani, R.; Jaubert, J.-N. Are safe results obtained when the PC-SAFT equation of state is applied to ordinary pure chemicals? Fluid Phase Equilib. 2010, 295, 76−92. (34) Alsaifi, N. M.; Al Aslani, I.; Binous, H.; Wang, Z.-G. A priori determination of the region of the three physical volume root loci in the Perturbed-Chain SAFT EOS. Fluid Phase Equilib. 2017, 434, 152−166. (35) Aslam, N.; Sunol, A. K. Reliable computation of all the density roots of the statistical associating fluid theory equation of state through global fixed-point homotopy. Ind. Eng. Chem. Res. 2006, 45, 3303−3310. (36) 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, 1718−1722. (37) Yelash, L.; Müller, M.; Paul, W.; Binder, K. A global investigation of phase equilibria using the perturbed-chain statistical-associating-fluid-theory approach. J. Chem. Phys. 2005, 123, 014908. (38) Yelash, L.; Müller, M.; Paul, W.; Binder, K. Artificial multiple criticality and phase equilibria: an investigation of the PC-SAFT approach. Phys. Chem. Chem. Phys. 2005, 7, 3728−3732. (39) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (40) Wertheim, M. S. Fluids with highly directional attractive forces. I. Statistical thermodynamics. J. Stat. Phys. 1984, 35, 19−34. (41) Wertheim, M. S. Fluids with highly directional attractive forces. II. Thermodynamic perturbation theory and integral equations. J. Stat. Phys. 1984, 35, 35−47. (42) Wertheim, M. S. Fluids with highly directional attractive forces. III. Multiple attraction sites. J. Stat. Phys. 1986, 42, 459−476. (43) Wertheim, M. S. Fluids with highly directional attractive forces. IV. Equilibrium polymerization. J. Stat. Phys. 1986, 42, 477−492. (44) Jackson, G.; Chapman, W. G.; Gubbins, K. E. Phase equilibria of associating fluids: Spherical molecules with multiple bonding sites. Mol. Phys. 1988, 65, 1−31. (45) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation-of-state solution model for associating fluids. Fluid Phase Equilib. 1989, 52, 31−38. (46) Ghonasgi, D.; Chapman, W. G. A new equation of state for hard chain molecules. J. Chem. Phys. 1994, 100, 6633−6639. (47) Alsaifi, N. M.; Patey, G. N.; Englezos, P. A general treatment of polar-polarizable systems for an equation of state. Chem. Eng. Res. Des. 2014, 92, 2936−2946.

(48) Al-Saifi, N. M.; Hamad, E. Z.; Englezos, P. Prediction of vapor−liquid equilibrium in water−alcohol−hydrocarbon systems with the dipolar perturbed-chain SAFT equation of state. Fluid Phase Equilib. 2008, 271, 82−93. (49) Zhang, P.; Alsaifi, N. M.; Wu, J.; Wang, Z.-G. Salting-Out and Salting-In of Polyelectrolyte Solutions: A Liquid-State Theory Study. Macromolecules 2016, 49, 9720−9730. (50) Zhang, P.; Alsaifi, N. M.; Wu, J.; Wang, Z.-G. Polyelectrolyte complex coacervation: Effects of concentration asymmetry. J. Chem. Phys. 2018, 149, 163303. (51) Dufal, S.; Lafitte, T.; Galindo, A.; Jackson, G.; Haslam, A. J. Developing intermolecular-potential models for use with the SAFTVR Mie equation of state. AIChE J. 2015, 61, 2891−2912. (52) Blas, F. J.; Vega, L. F. Prediction of Binary and Ternary Diagrams Using the Statistical Associating Fluid Theory (SAFT) Equation of State. Ind. Eng. Chem. Res. 1998, 37, 660−674. (53) Blas, F. J.; Vega, L. F. Thermodynamic behaviour of homonuclear and heteronuclear Lennard-Jones chains with association sites from simulation and theory. Mol. Phys. 1997, 92, 135−150. (54) Huang, S. H.; Radosz, M. Equation of state for small, large, polydisperse, and associating molecules. Ind. Eng. Chem. Res. 1990, 29, 2284−2294. (55) Chen, S. S.; Kreglewski, A. Applications of the Augmented van der Waals Theory of Fluids.: I. Pure Fluids. Berichte Bunsenges. Für Phys. Chem. 1977, 81, 1048−1052. (56) Binous, H.; Shaikh, A. A. Introduction of the arc-length continuation technique in the chemical engineering graduate program at KFUPM. Comput. Appl. Eng. Educ. 2015, 23, 344−351. (57) Al-Mubaiyedh, U. A.; Binous, H. Teaching arc-length continuation in the chemical engineering graduate program using MATLAB. Comput. Appl. Eng. Educ. 2018, 26, 1033. (58) Gallo, P.; et al. Water: A tale of two liquids. Chem. Rev. 2016, 116, 7463−7500. (59) Waseem, M. S.; Alsaifi, N. M. Prediction of vapor-liquidhydrate equilibrium conditions for single and mixed guest hydrates with the SAFT-VR Mie EOS. J. Chem. Thermodyn. 2018, 117, 223− 235. (60) Gubbins, K. E. Perturbation theories of the thermodynamics of polar and associating liquids: A historical perspective. Fluid Phase Equilib. 2016, 416, 3−17. (61) Colina, C. M.; Turrens, L. F.; Gubbins, K. E.; Olivera-Fuentes, C.; Vega, L. F. Predictions of the Joule−Thomson Inversion Curve for the n-Alkane Series and Carbon Dioxide from the Soft-SAFT Equation of State. Ind. Eng. Chem. Res. 2002, 41, 1069−1075. (62) Vega, L. F.; Llovell, F.; Blas, F. J. Capturing the Solubility Minima of n-Alkanes in Water by Soft-SAFT. J. Phys. Chem. B 2009, 113, 7621−7630. (63) Pàmies, J. C.; Vega, L. F. Vapor- liquid equilibria and critical behavior of heavy n-alkanes using transferable parameters from the soft-SAFT equation of state. Ind. Eng. Chem. Res. 2001, 40, 2532− 2543. (64) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The LennardJones equation of state revisited. Mol. Phys. 1993, 78, 591−618. (65) Aursand, P.; Gjennestad, M. A.; Aursand, E.; Hammer, M.; Wilhelmsen, Ø. The spinodal of single-and multi-component fluids and its role in the development of modern equations of state. Fluid Phase Equilib. 2017, 436, 98−112. (66) Kunz, O.; Wagner, W. The GERG-2008 wide-range equation of state for natural gases and other mixtures: an expansion of GERG2004. J. Chem. Eng. Data 2012, 57, 3032−3091. (67) Huang, S. H.; Radosz, M. Equation of state for small, large, polydisperse, and associating molecules: extension to fluid mixtures. Ind. Eng. Chem. Res. 1991, 30, 1994−2005. (68) Barker, J. A.; Henderson, D. Perturbation theory and equation of state for fluids. II. J. Chem. Phys. 1967, 47, 4714−4721. (69) Cotterman, R. L.; Schwarz, B. J.; Prausnitz, J. M. Molecular thermodynamics for fluids at low and high densities. Part I: Pure fluids containing small or large molecules. AIChE J. 1986, 32, 1787− 1798. M

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (70) Elliott, J. R., Jr; Suresh, S. J.; Donohue, M. D. A simple equation of state for non-spherical and associating molecules. Ind. Eng. Chem. Res. 1990, 29, 1476−1485. (71) Elliott, J. R.; Lira, C. T. Introductory Chemical Engineering Thermodynamics; Prentice Hall PTR: Upper Saddle River, NJ, 2012); p 802. (72) To generate Figure 12, we estimated the adjustable parameters for ethane using the same methodology applied in CK-SAFT.54 The adjustable parameters were found to be m = 1.5858, σ = 3.2469 Å, ε/k = 230.60 K by fitting liquid density and vapor pressure. For the effective hard sphere diameter, the temperature dependence is given by d = σ[1 − c exp(u/kT)], where c = 0.41052 and u/k = ε/k(1 − 1.189/T).

N

DOI: 10.1021/acs.iecr.8b04656 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX