# Clapeyron.jl: An extensible, open-source fluid-thermodynamics toolkit

Pierre J. Walker  Equally contributing authors<EMAIL_ADDRESS>Hon-Wa
Yew  Equally contributing authors<EMAIL_ADDRESS>Andrés Riedemann  [

###### Abstract

Thermodynamic models are often vital when characterising complex systems,
particularly natural gas, electrolyte, polymer, pharmaceutical and biological
systems. However, their implementations have historically been abstruse and
cumbersome, and as such, the only options available were black-box commercial
tools. In this article, we present Clapeyron.jl: a pioneering attempt at an
open-source fluid-thermodynamics toolkit to build and make use of
thermodynamic models. This toolkit is built in Julia, a modern language for
scientific computing known for its ease of use, extensibility, and first-class
support for differentiable programming. We currently support more models than
any package available, including standard cubic (SRK, PR, PSRK, etc.),
activity-coefficient (NRTL, UNIFAC, etc.), COSMO-based, and the venerable SAFT
equations. The property-estimation methods supported are extensive, including
bulk, VLE, LLE, VLLE and critical properties. With Clapeyron.jl, researchers
and enthusiasts alike will be able to focus on the application and worry less
about the implementation.

###### keywords:

Equations of State, Cubics, SAFT, Activity, COSMO-SAC, Open-source,
Thermodynamics, Julia

[California Institute of Technology] Division of Chemistry and Chemical
Engineering, California Institute of Technology, Pasadena, CA, 91125, U.S.A
[ICL]Department of Chemical Engineering, Imperial College London, SW7 2AZ,
United Kingdom Concepcion] Department of Chemical Engineering, University of
Concepción, Concepción, Región del Bio Bio, Chile Equally contributing authors

##  1 Introduction

Thermodynamic models, such as equations of state and activity-coefficient
models, are used to predict the properties of bulk matter at specified
conditions. Today, they have found applications in natural gas1, electrolyte2,
3, polymer4, 5, 6, metallic7, pharmaceutical8, 9 and biological10 systems.
They have become ubiquitous in modelling for both academia and industry1 and,
moreover, as the needs and demands placed upon this modelling grow ever more
larger, these models have become more sophisticated. The thermodynamic
modelling scene has now become quite complex, and it has become very difficult
to keep track of all developments in the various fields.

Figure 1: Condensed family tree of the thermodynamic models developed in the
past two centuries.

The family tree of thermodynamic models, illustrated in figure 1, begins in
1834, when Clapeyron 11 combined the laws of Boyle, Charles, Avogadro and Gay-
Lussac, resulting in what we now refer to as the ideal gas law, which is
arguably the most famous equation of state; indeed, it has become common
knowledge to almost anyone with scientific literacy. It would not be until
decades later that van der Waals 12 then introduced his equation of
state—which paved the way to an entire class of equations of state known as
the cubics. This class of equations proved to be extremely popular, primarily
due to their simple functional form. Some of the notable cubic equations of
state: Redlich–Kwong13 (RK), Soave–Redlich–Kwong14 (SRK) and Peng–Robinson15
(PR), continue to be the most-used in industry today. Subsequent modifications
of such equations have also been developed in the hopes of improving their
accuracy whilst maintaining their simple functional form.16, 17, 18, 19, 20,
21, 22, 23 These modifications alone are numerous and, for those unfamiliar
with the field, can be difficult to navigate.

The van der Waals’ equation of state also turned out to be the foundation of
another branch of models known as perturbation models, which are derived from
first principles using statistical mechanics. It was not until the development
of the Statistical Associating Fluid Theory (SAFT), introduced by Chapman et
al. 24, that such equations became a viable alternative. SAFT-type equations
retain a much-stronger physical basis than cubic equations, and include terms
to explicitly account for dispersive and associative interactions. The use of
SAFT-type equations is continually becoming more widespread; indeed, the
popularity of the PC-SAFT equation, developed by Gross and Sadowski 25, rivals
that of the cubic equations of state. This increasing uptake of SAFT-type
equations is a reflection of their broader applicability and more-predictive
nature – but this flexibility has come at the cost of greater complexity in
the equations which, together with their rapid development in recent decades,
has meant that they have become very difficult to access for people outside
the field.

In parallel to this tree, other types of thermodynamic models have also been
developed. These include multi-parameter equations of state such as IAPWS-9726
and GERG-200827 which are restricted to specific systems but with levels of
accuracy much higher than standard equations of state (we refer to them as
Empirical Helmholtz equations of state). Another branch is the activity-
coefficient models (such as the Wilson equation28, NRTL29 and UNIFAC30) which,
although not as flexible as typical equations of state, have become very
popular due to their simple form, large parameter databases and ease of
computing typically complicated properties. Whilst also technically part of
this branch, COSMO-based models31, such as COSMO-RS32 and COSMO-SAC33 have
proven to be very powerful, fully predictive thermodynamic models. Activity-
coefficient models have also been combined with cubic equations of state17, 21
to improve the modelling of mixtures. This combination has resulted in some of
the most-reliable equations of state, such as PSRK34 and VTPR35.

With such a wide selection of thermodynamic models, it is very easy to become
overwhelmed. Unfortunately, the difficulty with these models does not end
here. Typically, thermodynamic models provide a single property (usually the
pressure, Helmholtz free energy or activity coefficients) at specified
conditions (usually volume, temperature and composition). Obtaining other
useful properties can be an arduous process. For example, just trying to
obtain bulk properties such as the heat capacity will usually involve taking
multiple derivatives of sometimes very complicated functions. As another
example, obtaining vapour–liquid equilibrium (VLE) or critical properties will
usually involve solving non-linear systems of equations made up of these high-
order derivatives. Finally, flash routines, which solve for the phase
equilibrium at specified conditions, making them one of the most important
tools when using thermodynamic models, involve complicated, atypical
optimisation problems; indeed, such is the complexity of the problem that an
entire book 36 has been written on the subject.

With all these complicated aspects when using thermodynamic models, it
shouldn’t be surprising that, until recently, only commercial packages such as
ASPEN37 and gPROMS38 have provided the necessary tools to work with these
models. The difficulty with such packages, aside from the financial barrier,
is their ‘black-box’ nature where, having access only to the front-end, users
cannot modify the underlying code. Users are also limited to the selection of
models provided by the package.

However, in recent years, the push for open-source software has led to the
creation of many open-source thermodynamic modelling packages. Examples of
this include CoolProp39 and thermo40, although these place less emphasis on
equations of state and activity-coefficient models. For packages focusing
specifically on the latter, very recent examples include phasepy41 written in
Python, thermopack42 written in FORTRAN, FeOs43 written in RUST and teqp44
written in C++. We note also that, famously, the PC-SAFT equation was released
open-source when it was initially published, something that has unfortunately
not been repeated for subsequent SAFT equations.

Unfortunately, despite giving users complete access to source code, these
endeavours still suffer from a few problems. Firstly, like their commercial
counter-parts, they do not provide a wide range of thermodynamic models. Of
course, it should be possible for users to use the open-source code as a
template to implement their own models. However, as mentioned previously, in
order to use these models in a productive way, one requires access to higher-
order derivatives which will need to be obtained analytically in the case of
phasepy and thermopack. To see where the difficulty arises here, it is
recommended that readers look at thermopack’s documentation, which contains
these derivatives. Particularly for SAFT-type models, these derivatives become
very tedious.

In the case of FeOs and teqp, tools such as hyper dual numbers45 and automatic
differentiation46 are leveraged to avoid having to obtain derivatives
analytically. The difficulty with these packages is linked to the programming
language used; low-level languages such as Rust and C++ are chosen to maximise
performance, but they are less user friendly. To make it accessible, an
interpreted front-end (such as Python) is generally used, which helps with the
adoption, but makes it difficult for a typical user to make extensions to the
framework without substantial effort.

In this article we introduce Clapeyron.jl, a package named in honour of the
progenitor of equations of state, that has been designed to avoid these
difficulties. Like FeOs and teqp, Clapeyron.jl makes use of automatic
differentiation but, unlike these packages, is written entirely in the novel,
easy-to-use scientific-computing language, Julia. It already incorporates over
30 thermodynamic models, including all the well-known SAFT-type equations,
many cubic equations and, in addition, activity-coefficient models and COSMO-
SAC. Moreover, Clapeyron.jl is presented not just as open-source, but also as
a user-extensible package, with users encouraged to submit their own tools and
models for incorporation to the package, so that Clapeyron.jl can grow and
become more flexible as the field develops.

Note that all figures in this manuscript can be re-created using the jupyter
notebook file provided in the supplementary information. Whilst these figures
contain multiple thermodynamic models, it is not the intention of this article
to compare the performance of these models, only to illustrate specific
functionalities of Clapeyron.jl.

The remainder of the article is set out as follows: in section 2, we provide
details on Clapeyron.jl’s architecture, particularly model objects and
functions that are the center of all computations. In section 3, we given a
comprehensive overview of the thermodynamic models available in Clapeyron.jl
and some of the special features available in each of them. Subsequently, in
section 4, we describe how the various property estimation methods are
implemented in an extensible way. We then highlight some current developments
and future plans for the package in section 5. Finally, in section 6, we
summarise the features of Clapeyron.jl.

##  2 Clapeyron.jl Architecture

Figure 2: Architecture of the Clapeyron.jl package.

Clapeyron.jl is divided into four main components, illustrated visually in
figure 2:

  * •

database \- the logic for extraction of parameters from comma-separated value
files for the implementation into the model objects. More information on the
structure of our databases can be fund in the supplementary information.

  * •

models \- functions that produce the model object that is at the center of
Clapeyron.jl.

  * •

methods \- the thermodynamic methods to obtain thermodynamic properties from
the models.

  * •

solvers \- custom solvers provided by Clapeyron.jl.

This structure is designed to be easy to follow for brand-new users of the
package. One can examine each aspect of Clapeyron.jl on the GitHub repository
(see section 7). However, to appreciate how Clapeyron.jl works, one needs
understand two of its key aspects: the model object and the eos function.
These are discussed in detail in the subsequent sections.

Figure 3: An example of the abstract type hierarchy in the Clapeyron.jl.

###  2.1 Model object

In Clapeyron.jl, thermodynamic models are created as structs (a composite
type), which contain all of the parameters needed to completely identify the
species and derive the thermodynamic properties. These structs are created as
sub-types of an abstract type in order to create a hierarchy in which the
model objects reside. Functions that are related to the models are written in
the corresponding Julia script that dispatch on instances of these abstract
types.

As an example, if we wanted to use the PC-SAFT equation of state for a water +
carbon dioxide system, we would construct the model object by calling the
constructor as:

⬇

julia> model = PCSAFT([”water”,”carbon dioxide”])

PCSAFT{BasicIdeal} with 2 components:

”water”

”carbon dioxide”

Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol

\parjulia> model.params.epsilon

PairParam{PairParam{Float64}}[”water”, ”carbon dioxide”]) with values:

2×2 Matrix{Float64}:

366.51 250.508

250.508 169.21

\parjulia> model.references

2-element Vector{String}:

”10.1021/ie0003887”

”10.1021/ie010954d”

As we can see above, the model object contains all the parameters that
characterise the water+carbon dioxide system. We also include the references
to the source of the parameters and thermodynamic model, ensuring that the
developers of these equations receive the appropriate credit.

This model also has a type associated with it. In Clapeyron.jl, we have
developed a hierarchy of types (a condensed version is illustrated in figure
3), the head of which is the EoSModel abstract type. From this parent type, we
divide the hierarchy into SAFTModel, CubicModel, ActivityModel, and so forth.
The SAFTModel type is then further divided into types associated with each
SAFT-type equation of state. This model object is a sub-type of the
PCSAFTModel type. What this effectively means is that Julia will only perform
operations on model using the functions whose input most closely resembles its
type.

To illustrate this we consider another model object constructed from the
sPCSAFT function (i.e. model2=sPCSAFT(["water","carbon dioxide"])) and that,
accordingly, will be associated with the sPCSAFTModel type (which is a sub-
type of PCSAFTModel). Only three equations are different between the theories
expressed in the PC-SAFT and sPC-SAFT equations of state. One of these is the
reduced Helmholtz free energy for a hard-sphere fluid. Within Clapeyron.jl,
the function corresponding to this equation is a_hs, written for each equation
of state as:

⬇

\par# PC-SAFT code

function a_hs(model::PCSAFTModel, V, T, z,_data=@f(data))

_,ζ0,ζ1,ζ2,ζ3,_ = _data

return 1/ζ0 * (3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3))

end

\par# sPC-SAFT code

function a_hs(model::sPCSAFTModel, V, T, z,_data=@f(data))

_,_,_,_,η,_ = _data

return (4η-3η^2)/(1-η)^2

end

Calling a_hs using the model object will tell Julia to only use the first
function, whilst calling a_hs using the model2 will tell Julia to use the
latter function. This powerful feature allows developers to define variants of
thermodynamic models (like sPC-SAFT) with very few lines of code, and even
define unique methods for specific models, should the need arise.

The above is also an example of how Clapeyron.jl makes full use of the Unicode
characters, bringing the code much closer to how the corresponding equations
are presented in their original publications. The use of Unicode characters is
strongly supported within the Julia community.

####  2.1.1 Group contribution

Clapeyron.jl also supports group-contribution (GC) models, where parameters
are associated with groups which are combined to represent species. Some GC
models that we currently support include SAFT-$\gamma$ Mie47 equation of
state, the UNIFAC48 activity-coefficient model and, Walker’s model49 and the
Joback Method50 for ideal-gas heat capacities. In order to form species from
groups within Clapeyron.jl, one can specify the type and number of groups
within a species as using a tuple containing the name of the species, and a
list of key-value pairs where the key is the name of the group and the value
is the group multiplicity:

⬇

julia> model = SAFTgammaMie([(”ethanol”,[”CH3”=>3,”CH2OH”=>1]),(”ethyl
acetate”, [”CH3”=>2,”COO”=>1,”CH2”=>1])])

SAFTgammaMie{BasicIdeal} with 2 components:

”ethanol”: ”CH3”=>3, ”CH2OH”=>1

”ethyl acetate”: ”CH3”=>2, ”COO”=>1, ”CH2”=>1

Contains parameters: segment, mixedsegment, shapefactor, lambda_a, lambda_r,
sigma, epsilon, epsilon_assoc, bondvol

This is in a similar vein to how species are constructed from groups in tools
such as gPROMS38.

###  2.2 eos function

With a large number of thermodynamic models available within a single package,
it is necessary to ensure that all models have a common structure at some
level. In Clapeyron.jl, we dispatch on a function called eos, which is the
Helmholtz free energy. This function is at the center of all computations in
Clapeyron.jl. For most models, we would generally make use of the built-in
function that dispatches on the top-level EoSModel, where it is made up of an
ideal and residual contribution (since most thermodynamic models are built
this way):

⬇

# General function

function eos(model::EoSModel, V, T, z=SA[1.0])

return N_A*k_B*Σ(z)*T * (a_ideal(idealmodel(model),V,T,z)

+a_res(model,V,T,z))

end

\par# Ideal models only

function eos(model::IdealModel, V, T, z=SA[1.0])

return N_A*k_B*Σ(z)*T *a_ideal(model,V,T,z)

end

Here V is the total volume (in m3), T is the temperature (in K) and z is the
amount of components (in mol).

It is, of course, possible to express the Helmholtz free energy in a different
way. For example, the if model is a purely ideal model, there is no residual
contribution. Using multiple dispatch, we can define an eos function specific
to ideal models (using the IdealModel type).

##  3 Models

Table 1: List of supported thermodynamic models within Clapeyron.jl, along
with their function names. Models that are sub-types of a main model that are
also supported are indicated as variants.

Type | Main | Function | Variant | Function
---|---|---|---|---
Cubics | van der Waals12 | vdW |  |
Redlich-Kwong13 | RK | Soave-Redlich-Kwong14 | SRK
Predictive SRK51 | PSRK
Peng-Robinson15 | PR | PR-7852 | PR78
UMR-PR53 | UMRPR
Volume-Translated PR35 | VTPR
| SAFT54, 24 | ogSAFT |  |
| CK-SAFT55, 56 | CKSAFT | simplified CK-SAFT57 | sCKSAFT
| BACK-SAFT58 | BACKSAFT |  |
| LJ-SAFT59, 60 | LJSAFT |  |
| soft-SAFT61 | softSAFT |  |
| CPA62, 63 | CPA | simplified CPA64 | sCPA
| SAFT-VR SW65 | SAFTVRSW |  |
| PC-SAFT25, 4 | PCSAFT | simplified PC-SAFT66 | sPCSAFT
| SAFT-VR Mie67, 68 | SAFTVRMie | SAFT-VRQ Mie69 | SAFTVRQMie
SAFT | SAFT-$\gamma$ Mie47, 70 | SAFTgammaMie |  |
Activity | Wilson28 | Wilson |  |
NRTL29 | NRTL |  |
UNIQUAC71 | UNIQUAC |  |
UNIFAC | N/A | og-UNIFAC30 | ogUNIFAC
mod-UNIFAC48 | UNIFAC
COSMO-SAC | N/A | COSMO-SAC-0233 | COSMOSAC02
COSMO-SAC-1072 | COSMOSAC10
COSMO-SAC-dsp73 | COSMOSACdsp
| Reid Polynomial74 | ReidIdeal |  |
| Joback Method50 | JobackIdeal |  |
Ideal | Walker Model49 | WalkerIdeal |  |
Empirical Helmholtz | IAPWS-9526 | IAPWS95 |  |
GERG-200827 | GERG2008 |  |
Propane Reference75 | PropaneRef |  |
SPUNG | SPUNG76 | SPUNG |  |

Clapeyron.jl supports over 30 benchmarked thermodynamic models, most of which
are listed in table 1. The highly flexible nature of Clapeyron.jl enables us
to express many thermodynamic models in a modular way, which allows users to
mix-and-match different components in an unprecedented manner.

###  3.1 Cubic equations of state

Table 2: List of supported $\alpha$-functions, volume translation and mixing rule methods wiht Clapeyron.jl, along with their function names. Component | Name | Function
---|---|---
Alpha function | Redlich–Kwong13 | RKAlpha
Soave14 | SoaveAlpha
Peng–Robinson15 | PRAlpha
PR-7852 | PR78Alpha
Boston-Matthias18 | BMAlpha
Twu22 | TwuAlpha
Magoulas-Tassios77 | MTAlpha
Volume Translation | Rackett16 | RackettTranslation
Peneloux20 | PenelouxTranslation
Magoulas-Tassios77 | MTTranslation
Mixing rule | van der Waals one-fluid78 | vdW1fRule
Kay79 | KayRule
Huron-Vidal17 | HVRule
Modified Huron-Vidal first order21 | MHV1Rule
Modified Huron-Vidal second order21 | MHV2Rule
LCVM80 | LCVMRule
Wong-Sandler81 | WSRule

In general, cubics are written in the following form:

$p=\frac{Nk_{\text{B}}T}{V-Nb}+\frac{a(T)}{(V-Nb_{1})(V-Nb_{2})}\,,$ (1)

where $p$, $N$, $V$, $T$ and $k_{\text{B}}$ are the pressure, total number of
particles, volume, temperature and the Boltzmann constant, respectively.
$a(T)$ and $b$ are related to the inter-molecular interactions and excluded
volume of the system, respectively. $b_{1}$ and $b_{2}$ are related to the
parameter $b$ by a certain factor dependent on the equation of state.

(a) Saturation curve of water for different $\alpha$ functions and volume
translations.82

(b) VLE envelope of water for different $\alpha$ functions and volume
translations.82

(c) $pxy$ diagram of CO (1) + CO2 (2) at 218.15 K for different $\alpha$
functions and mixing rules.83

(d) $pxy$ diagram of methanol (1) + benzene (2) at 433.15 K for different
mixing rules.84

Figure 4: Examples of the different capabilities of Clapeyron.jl for cubic
equations of state.

In contrast to many of the other models listed in table 1, cubic equations of
state offer a great deal of flexibility in three particular aspects (these
have been summarised in table 2):

  * •

$\alpha$-functions: Initially introduced by Soave 14 to improve the
predictions made by the Redlich–Kwong equation13, numerous $\alpha$-functions
have been developed to improve the performance of the equations of state. The
$\alpha$-functions typically appear within the function $a(T)$ shown in
equation 1 as:

$a(T)=a\alpha(T)\,,$ (2)

where $a$ is a parameter usually related to the critical point of a species
and $\alpha(T)$ is the $\alpha$-function. To give an example of this, the
$\alpha$-function within the Redlich–Kwong equation is:

$\alpha(T)=\frac{1}{\sqrt{T_{r}}}\,,$ (3)

where $T_{r}$ is the reduced temperature relative to the critical point. The
Soave–Redlich–Kwong equation differs from the RK equation through its
$\alpha$-function:

$\alpha(T)=\left(1+(0.480+1.547\omega-0.176\omega^{2})\left(1-\sqrt{T_{r}}\right)\right)^{2}\,,$
(4)

where $\omega$ is the acentric factor. Within Clapeyron.jl, it is possible to
switch between $\alpha$-functions by simply providing an additional optional
argument. This is, in fact, how we construct the Soave–Redlich–Kwong equation
(rather than writing an entirely separate section of code):

⬇

julia> model = SRK([”butane”])

RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule} with 1 component:

”butane”

Contains parameters: a, b, Tc, Pc, Mw

\parjulia> model = SRK([”butane”];alpha=BMAlpha)

RK{BasicIdeal, BMAlpha, NoTranslation, vdW1fRule} with 1 component:

”neon”

Contains parameters: a, b, Tc, Pc, Mw

It is equally as easy to switch to another $\alpha$-function listed in table
2, as shown above. Examples of this are shown in figures 4(a) and 4(b), where
switching $\alpha$-functions can improve the accuracy of our results.
Furthermore, if one needs to model a species beyond its critical point, the
Boston–Mathias $\alpha$-function can be specified, which can also have an
impact on VLE properties as shown in figure 4(c).

  * •

Volume translation: One often finds that the cubic equations perform poorly
when predicting liquid densities. To improve upon this without affecting the
simple functional form of the equation, one can ‘translate’ the volume used in
the equation of state by a certain value $c$:

$V=V_{\text{EoS}}+Nc\,,$ (5)

where $V_{\text{EoS}}$ represents the untranslated volume in the equation of
state. There are a few different methods which are listed in table 2 and, like
the $\alpha$-functions, it is very easy to switch between them for a
particular cubic equation:

⬇

julia> model = PR([”water”])

PR{BasicIdeal, PRAlpha, NoTranslation, vdW1fRule} with 1 component:

”water”

Contains parameters: a, b, Tc, Pc, Mw

julia> volume(model,1e5,298.15)

2.1390770714552146e-5

julia> model = PR([”water”];

translation=RackettTranslation)

PR{BasicIdeal, PRAlpha, RackettTranslation, vdW1fRule} with 1 component:

”water”

Contains parameters: a, b, Tc, Pc, Mw

julia> volume(model,1e5,298.15)

1.8705398679867126e-5

We illustrate the capability of volume translation in Clapeyron.jl in figures
4(a) and 4(b); as shown in these figures, whilst it doesn’t have an impact on
the saturation curve, it can improve the predicted liquid densities.

  * •

Mixing rules: When one wishes to model mixtures with a cubic equation of
state, mixing rules are required to obtain the parameters $a$ and $b$ in
equation 1. The most famous and simplest of these is the van der Waals one-
fluid mixing rule:

$\bar{a}=\sum_{i}\sum_{j}x_{i}x_{j}a_{ij}\,,$ (6)
$\bar{b}=\sum_{i}\sum_{j}x_{i}x_{j}b_{ij}\,,$ (7)

where $x_{i}$ denotes the composition of species $i$, the over bars indicate
the parameters characterising the one-fluid average and subscripts $ij$ denote
parameters characterising the interactions between species $i$ and $j$. When
$i=j$ we refer to the familiar species-specific parameters; when $i\neq j$, we
are referring to the unlike parameters which are typically obtained from
combining rules:

$a_{ij}=\sqrt{a_{ii}a_{jj}}(1-k_{ij})\,,$ (8)
$b_{ij}=\frac{b_{ii}+b_{jj}}{2}\,,$ (9)

where $k_{ij}$ is the binary interaction parameter between $i$ and $j$. The
van der Waals one-fluid mixing rule provides satisfactory results when dealing
with similar species. For more-complicated systems, advanced mixing rules
should be used. Examples of these are listed within table 2; details on the
formulation of these can be found in their respective references. What sets
these mixing rules apart from the van der Waals one-fluid (and Kay) rule is
that an activity coefficient model must be specified as well. The flexibility
of Clapeyron.jl allows one to choose between any of the activity-coefficient
models listed in table 1. For example:

⬇

julia> model = PR([”methanol”,”benzene”];mixing=KayRule)

PR{BasicIdeal, PRAlpha, NoTranslation, KayRule} with 2 components:

”methanol”

”benzene”

Contains parameters: a, b, Tc, Pc, Mw

julia> model = PR([”methanol”,”benzene”];mixing=HVRule, activity=UNIFAC)

PR{BasicIdeal, PRAlpha, NoTranslation, HVRule{UNIFAC}} with 2 components:

”methanol”

”benzene”

Contains parameters: a, b, Tc, Pc, Mw

One can observe, in figures 4(c) and 4(d), the impact of using different
mixing rules on the predicted VLE behaviour of a mixture.

As one can imagine, with the flexibility to switch between $\alpha$-functions,
volume-translation methods and mixing rules with nothing more than a change in
an input argument, Clapeyron.jl is able to provide a variety of cubic
equations of state with a relatively small amount of code. An example of the
benefits of this is the ability to ‘create’ some of the high-accuracy cubic
equations of state such as Predictive-SRK51 (example calculations illustrated
in figure 4(c)) and Volume-Translated-PR35 (exemplified in figure 4(d)).

###  3.2 SAFT equations of state

(a) VLE envelope of methanol for using different SAFT equations.82

(b) $pxy$ diagram of CO2 (1) + water (2) at 333.15 K using different SAFT
equations.85

Figure 5: Examples of the various SAFT equations available in Clapeyron.jl.

Clapeyron.jl provides the most extensive list of SAFT-type equations of state
in any open-source or commercial package available. As an example, most of the
supported SAFT equations are included in figure 5(a); methanol was selected as
a species here as it requires the use of all three terms in the SAFT equation.
For future reference, we remind the readers of the general formulation of the
SAFT equation24:

$A_{\text{res.}}=A_{\text{seg.}}+A_{\text{chain}}+A_{\text{assoc.}}\,,$ (10)

where $A$ is the Helmholtz free energy and the subscripts res., seg., chain
and assoc. refer to the residual contribution, the contribution due to the
presence of segments, the contribution to the formation of chains of segments
and the contribution due to association between segments, respectively. Whilst
most SAFT equations, at the surface, appear the same, this is where most of
the similarities end. The one commonality between all SAFT equations is the
association term. As a result, we have common association solver used by all
SAFT-type models. Here, the association fraction of sites $a$ on species $i$,
$X_{ia}$, is solved for using the following mass-action equation:

$X_{ia}=\frac{1}{1+\sum_{j,b}{n_{jb}X_{jb}\Delta_{ijab}}\rho z_{i}}\,,$ (11)

where $\Delta_{ijab}$ is the association strength between site $a$ in
component $i$ and site $b$ in component $j$. As this expression is identical
in all SAFT equations, all of them can use one common solver. Our
implementation requires only the definition of the association strength. In
the general case, the mass-action equation is solved successively until
convergence. However, analytical solutions do exist for special cases, for
example, when there is only one pair of sites.

Aside from the common association term, all other terms in SAFT equations are
too distinct to offer the same level of flexibility present in the cubic
equations. Nevertheless, we are able to provide certain levels of flexibility
to some of the SAFT equations:

  1. 1.

Variants: This feature is not exclusive to SAFT-type equations and we have
given an example in section 2. However, it has been used to provide variants
of SAFT-type equations for multiple base SAFT equations. These have been given
in table 1. There is one other equation of state we treat as a ‘variant’,
although not in the same sense as the example given previously. It is possible
to treat SAFT-$\gamma$ Mie as a variant of SAFT-VR Mie as they share
practically the same equations67, 47 only adapted for a group-contribution
method. To take advantage of this fact, within our implementation of
SAFT-$\gamma$ Mie, we map this equation to the SAFT-VR Mie equation in the
following way (we refer to the mapped model as vrmodel):

     * •

Segment term: Within the segment term of SAFT-$\gamma$ Mie, groups are treated
as separate entities. As a result one can simply treat the groups as species
when mapping from SAFT-$\gamma$ Mie to SAFT-VR Mie.

     * •

Chain term: Since the groups are combined to make a single chain in
SAFT-$\gamma$ Mie, the group-based parameters are averaged to obtain effective
species-based parameters. Aside from this averaging, the formulation of the
chain terms in both equations is identical. As such, one can simply feed the
effective parameters straight into the SAFT-VR Mie chain term.

     * •

Association term: This mapping is a little more difficult as variables are
related to the species, groups and sites. To handle this, the groups are
effectively ignored and the sites that are present are then linked to the
species. This means that we must introduce new labels for the association
sites, rather than just the standard ‘hydrogen’ (H) or ‘electron’ (e1 or e2)
sites67, 47; these labels now contain the groups on which the sites are
present. As an example:

⬇

julia> model = SAFTgammaMie([

(”ethanol”,[”CH3”=>3,”CH2OH”=>1]),

(”glycolic acid”, [”COOH”=>1,”CH2OH”=>1])]);

\parjulia> model.sites

SiteParam with 3 components:

”CH3”: (no sites)

”CH2OH”: ”H” => 1, ”e1” => 2

”COOH”: ”e2” => 2, ”H” => 1, ”e1” => 2

\parjulia> model.vrmodel.sites

SiteParam with 2 components:

”ethanol”: ”CH2OH{H}” => 1,

”CH2OH{e1}” => 2

”glycolic acid”: ”CH2OH{H}” => 1,

”CH2OH{e1}” => 2,

”COOH{e2}” => 2,

”COOH{H}” => 1,

”COOH{e1}” => 2

We can see that three groups with three site types have been mapped to two
species with five site types, making the model compatible with the SAFT-VR Mie
equation.

This unconventional use of multiple-dispatch has allowed for a simpler
implementation of what is usually accepted as one of the most-complex SAFT
equations.

  2. 2.

Cubic-plus-association (CPA): The CPA equation of state62 is rather unique, in
contrast to other SAFT-type equations. Whilst it retains the association term,
the segment and chain terms are replaced by the residual Helmholtz free energy
obtained from a cubic equation of state (in its original publication, SRK is
used to provide this contribution). Due to the flexibility of the cubic
equations of state, in literature, this flexibility is often extended to the
CPA equation. Similarly, within Clapeyron.jl, we have extended the flexibility
available for cubics to the CPA equation of state. For example:

julia¿ model = CPA([”water”,”ethanol”]) CPABasicIdeal, RKBasicIdeal, CPAAlpha,
NoTranslation, vdW1fRule with 2 components: ”water” ”ethanol” Contains
parameters: a, b, c1, Tc, epsilon_assoc, bondvol, Mw

julia¿ model = CPA([”water”,”ethanol”]; cubicmodel=PR,
translation=RackettTranslation, mixing=MHV1Rule, activity=UNIFAC)
CPABasicIdeal, PRBasicIdeal, CPAAlpha, RackettTranslation, MHV1RuleUNIFAC with
2 components: ”water” ”ethanol” Contains parameters: a, b, c1, Tc,
epsilon_assoc, bondvol, Mw This has resulted in possibly the most-flexible
implementation of the CPA equation of state, even including commercial
packages. For example, as illustrated in figure 5(b), one can switch the cubic
equation that is used in simplified CPA (a variant of CPA).

###  3.3 Activity-coefficient models

(a) $pxy$ diagram of ethanol (1) + water (2) at 433.15 K.86

(b) $Txy$ diagram of ethanol (1) + water (2) at 1.01 bar.87

Figure 6: Examples of the different capabilities of Clapeyron.jl for activity-
coefficient models.

Activity-coefficient models are a unique form of thermodynamic models, in
contrast to the equations of state discussed previously. Whilst we have
already shown them used within mixing rules to improve the predictions made by
cubic equations of state in multi-component mixtures, they can be used
directly to give us the activity coefficients of a mixture using which we can
perform vapour–liquid equilibrium calculations through a modified Raoult’s
law:

$y_{i}p=x_{i}\gamma_{i}p_{\text{sat},i}\,,$ (12)

where $p_{\text{sat},i}$ is the saturation pressure of component $i$.
Unfortunately, activity-coefficient models cannot obtain the saturation
pressure of individual components themselves; this must be provided by another
model. To enable this, all activity-coefficient models have an additional
input to select which equation of state will be used to determine the pure-
component saturation pressure:

⬇

julia> model = NRTL([”methanol”,”benzene”];

puremodel=PR)

NRTL{PR} with 2 components:

”methanol”

”benzene”

Contains parameters: a, b, c, Mw

Any equation of state can be used to obtain the saturation pressure, as long
as all component parameters are available; this is illustrated in figure 6(a).
Unfortunately, many of the limitations of activity-coefficient models remain:

     * •

Sub-critical regions only: Given activity-coefficient models rely on Raoult’s
law to obtain vapour–liquid equilibrium properties, it is required that all
components in the system be below their critical point. Thus, activity model
calculations are limited to the lowest critical temperature of the components
in the mixture.

     * •

Bulk properties: activity-coefficient models are solely designed for
vapour–liquid equilibria calculations and, technically, cannot be used to
obtain bulk properties. However, using Dalton’s law as an approximation, one
could obtain an expression for the total Helmholtz free energy of the system:

$A(V,T,\mathbf{z})=G^{E}(T,\mathbf{z})+G^{I}(T,\mathbf{z})+\sum_{i}z_{i}(G_{i}(V,T)-p_{i}(V,T)V)\,.$
(13)

where $G^{E}$, $G^{I}$, $G_{i}$ and $p_{i}$ are the excess Gibbs free energy,
ideal Gibbs free energy of mixing, Gibbs free energy of pure component $i$ and
pressure of pure component $i$, respectively. Although this provides a route
to the provision of bulk properties, it is, of course, a very crude
approximation and is not recommended.

Nevertheless, with the level of flexibility provided by Clapeyron.jl, one can
produce interesting comparisons between various activity-coefficient models
when they are used by themselves or within a mixing rule (as shown in figure
6(b)).

###  3.4 Ideal gas models

Often overlooked as they have no impact on equilibrium properties, the ideal
contribution to any equation of state can represent the dominant contribution
to various properties, particularly second-derivative properties.49 As such,
within Clapeyron.jl, we’ve provided both basic and more-advanced ideal
contribution models.

(a) Isobaric heat capacity of pentane at
5 MPa.82

(b) Joule–Thomson coefficient of CO2 at 5 MPa.75

Figure 7: Examples of different applications of the ideal term within
Clapeyron.jl.

Firstly, astute readers may have noticed that, whenever a model is created,
one of the sub-types it contains is BasicIdeal. This is the default ideal
model which requires no species-based parameters. It is given by the following
equation:

$A_{\text{ideal}}=Nk_{\text{B}}T\sum_{i}x_{i}\left[\log{\left(\frac{x_{i}}{VT^{3/2}}\right)}-1\right]\,.$
(14)

This equation is only meant to account for the translation modes of species,
thus, will give quantitatively inaccurate results for bulk properties such as
the heat capacities (see the results obtained from SRK in figure 7(a) and
SAFT-VR Mie in figure 7(b)), except for molecules (such as noble gases) where
only translation is important. To remedy this, we have provided three
alternative ideal terms: Reid polynomials74, the Joback method50 and Walker’s
model49. All three of these will account for the rotational and vibrational
modes missing from the BasicIdeal model. For example:

⬇

julia> model=SAFTVRMie([”benzene”])

SAFTVRMie{BasicIdeal} with 1 component:

”benzene”

Contains parameters: segment, sigma, lambda_a, lambda_r, epsilon,
epsilon_assoc, bondvol, Mw

\parjulia> isobaric_heat_capacity(model,1e5,298.15)

64.24460509411827

\parjulia> model=SAFTVRMie([”benzene”];idealmodel=WalkerIdeal)

SAFTVRMie{WalkerIdeal} with 1 component:

”benzene”

Contains parameters: segment, sigma, lambda_a, lambda_r, epsilon,
epsilon_assoc, bondvol, Mw

\parjulia> isobaric_heat_capacity(model,1e5,298.15)

125.6270421887885

Note that, in the above, we have combined a non-group-contribution residual
model with a group-contribution ideal model. Such combinations are allowed in
Clapeyron.jl as long as the groups for that particular species are available
within Clapeyron.jl. The opposite is also possible (combining a group-
contribution residual model with a non-group-contribution ideal model). Both
of these situations are illustrated in figures 7(a) and 7(b) where readers can
see that, even for properties like the Joule–Thomson coefficient, the ideal
model can have a significant impact on the overall property. Figure 7(a) also
highlights that one can obtain results from the ideal models alone in
Clapeyron.jl.

Whilst the current library of ideal terms is smaller than the other models,
the terms that are included provide extensive databases of parameters which
should be compatible with almost any species. If more ideal terms are needed,
these will be added over time or can be provided by the users themselves.

###  3.5 Empirical Helmholtz equations

What sets these equations apart from the rest is that they do not have a
common structure between them and are designed with particular systems in
mind. They are designed to be highly accurate for those particular systems
(see figures 7(a) and 9(a) as examples for GERG-200827).

The lack of a common structure and the correlative nature of these equations
prevents us from developing flexible frameworks in which to use them (save for
the SPUNG framework). Nevertheless, they serve as an example that
thermodynamic models within Clapeyron.jl need not belong to the cubic, SAFT or
activity-coefficient family of models. They can be stand-alone and still be
compatible with Clapeyron.jl’s solvers.

####  3.5.1 SPUNG

The State Research Program for Utilization of Natural Gas (SPUNG) framework
was initially devised by Jørstad 76 based on the work by Mollerup 88 using the
principle of corresponding states. The general idea here is, if we have a
model that is very good at representing a specific system, according to the
principle of corresponding states, we should be able to map the results of
this model to another similar model and expect to achieve similarly accurate
results. The details of this are explained in the original works. As one can
see in figure 7(a), if the two species are dissimilar, such as carbon dioxide
and propane, the results obtained can be quite poor.

This framework requires two things: 1. A highly accurate reference model for
the reference system; 2. A ‘shape’ model, which is used to map the reference
system to the desired system. Taking the following example:

⬇

julia> model = SPUNG([”propylene”])

Extended Corresponding States model

reference model: PropaneRef()

shape model: RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule}(”propylene”)

\parjulia> sat_pure(model,250)

(299496.50354575133, 7.191991655304882e-5, 0.006575782239022821)

As one can see, the default reference and shape models are the propane
reference and SRK equations of state. Where we begin seeing the true
flexibility of Clapeyron.jl is that the reference and shape equations can be
any of the models available, leading to an unimaginable number of
possibilities:

⬇

# SPUNG(components,refmodel,shapemodel,shaperef)

julia> model = SPUNG([”methanol”],

IAPWS95(),

PCSAFT([”methanol”]),

PCSAFT([”water”]))

Extended Corresponding States model

reference model: IAPWS95(”water”)

shape model: PCSAFT(”methanol”)

\parjulia> model = SPUNG([”water”],vdW([”methane”]),

SAFTgammaMie([”water”]),

SAFTgammaMie([”methane”]))

Extended Corresponding States model

reference model: vdW(”methane”)

shape model: SAFTgammaMie(”water”)

The latter example is bound to be inaccurate but helps illustrate just how
flexible Clapeyron.jl can be.

###  3.6 User-defined equations

The last class of thermodynamic models supported in Clapeyron.jl are those
defined by the users themselves. To limit the barrier of entry to novice
users, we have provided generic templates within our documentation. However,
this can summarised to just four components:

⬇

# Defining the parameters used by the model

struct NewEoSParam <: EoSParam

# Insert parameter names here along with their type (SingleParam, PairParam or
AssocParam)

end

\par# Defining an abstract type for this model type

abstract type NewModel <: SAFTModel end

\par# Creating a model struct called New, which is a sub-type of NewEoSModel,
and uses parameters defined in NewEoSParam

@newmodel NewEoS NewEoSModel NewEoSParam

\par# Function that will create the model object

function NewEoS(components; idealmodel=BasicIdeal, userlocations=String[],
ideal_userlocations=String[], verbose=false)

# Obtain a Dictionary of parameters. We pass in custom locations through the
optional parameter userlocations.

params = getparams(components; userlocations=userlocations, verbose=verbose)

\par# Insert any necessary transformation to the parameters (e.g. combining
rules)

\par# Create the parameter struct that we have defined.

packagedparams = NewEoSParam(”transformed parameters go here”)

\par# Although optional, it’s generally good practise to cite your models!

references = [”DOI 1”, ”DOI 2”]

\parmodel = NewEoS(packagedparams, idealmodel;
ideal_userlocations=ideal_userlocations, references=references,
verbose=verbose)

\par# Return the model object that you have just created.

return model

end

\par# Defining the equations

function a_res(model::NewEoSModel, V, T, z)

# Define all necessary equations that output the residual Helmholtz free
energy here

return a_res

end

However, this is assuming the desired model is entirely separate from all
those defined in Clapeyron.jl. If the user wishes to define a variant of an
existing thermodynamic model in Clapeyron.jl, the only equations they would
need to define are those that are different from the original equation coded
in Clapeyron.jl. This is a similar concept to how simplified PC-SAFT was
implemented in Clapeyron.jl, only needing to redefine the modified equations.
A full example of this is given in supplementary material for the SAFT-VR
Mie+AT equation developed by Walker et al. 89, requiring less than 150 lines
to implement (compared to the 900 lines needed to code the SAFT-VR Mie
equation).

In general, to assure compatibility with Clapeyron.jl, it is recommended to
simply modify existing equations in Clapeyron.jl or, at the very least, have
the new equation be a sub-type of CubicModel, SAFTModel or ActivityModel. This
increases the chance that the existing methods in Clapeyron.jl will instantly
be compatible with the new model. In the following section we illustrate why
this is the case.

##  4 Supported Methods

In this section we highlight some of the supported methods in Clapeyron.jl
which are listed in table 3. These are intended to cover a wide range of
properties whilst also being intuitive to use. Many of the standard properties
are provided although some unusual ones are available, such as the critical
points for both pure and mixture systems. These are rarely seen in commercial
packages such as gPROMS38 or ASPEN37. Details of their implementation will be
given in subsequent sections. However, all of these properties are compatible
with all models given in table 1 (with obvious exceptions such activity models
not being compatible with single-component VLE properties or critical
properties).

We note that the methods currently supported in Clapeyron.jl may not be as
robust as those provided in commercial packages, and we will point out some of
the short-comings of these methods, as well as some potential solutions to
these short-comings. However, as Clapeyron.jl is an open-source initiative, we
hope others in the field will be generous enough to contribute to the package,
continuously improving its methods.

All units are stated in MKS units by default. However, Clapeyron.jl is also
fully compatible with the Unitful.jl package90. We are effectively able to
inject units to the inputs, which Clapeyron.jl will automatically convert into
MKS units, and obtain the outputs with the correct units; this may be useful
when this package is used as part of a larger system with Unitful.jl support.
For example:

⬇

julia> using Unitful

\parjulia> import Unitful: bar, °C, mol

\parjulia> model = PCSAFT([”methanol”,”ethanol”])

PCSAFT{BasicIdeal} with 2 components:

”methanol”

”ethanol”

Contains parameters: Mw, segment, sigma, epsilon, epsilon_assoc, bondvol

\parjulia> Cp = isobaric_heat_capacity(model, 5bar, 25°C, [0.5mol,
0.5mol];output=cal*K^-1)

15.091222228910153 cal K⁻¹

Table 3: List of supported properties within Clapeyron.jl, along with their
function names.

Type | Property | Function
---|---|---
Bulk | Volume | V=volume(model, p, T, z)
Entropy | S=entropy(model, p, T, z)
Chemical Potential | μ=chemical_potential(model, p, T, z)
Internal Energy | U=internal_energy(model, p, T, z)
Enthalpy | H=enthalpy(model, p, T, z)
Gibbs Free Energy | G=gibbs_free_energy(model, p, T, z)
Helmholtz Free Energy | A=helmholtz_free_energy(model, p, T, z)
Isochoric Heat Capacity | Cv=isochoric_heat_capacity(model, p, T, z)
Isobaric Heat Capacity | Cp=isobaric_heat_capacity(model, p, T, z)
Isothermal Compressibility | βT=isothermal_compressibility(model, p, T, z)
Isentropic Compressibility | βS=isentropic_compressibility(model, p, T, z)
Speed of Sound | u=speed_of_sound(model, p, T, z)
Isobaric Expansivity | αV=isobaric_expansivity(model, p, T, z)
Joule–Thomson Coefficient | μJT=joule_thomson_coefficient(model, p, T, z)
Fugacity Coefficient | φ=fugacity_coefficient(model, p, T, z)
Activity Coefficient | γ=activity_coefficient(model, p, T, z)
Compressibility Factor | Z=compressibility_factor(model, p, T, z)
Inversion Temperature | T=inversion_temperature(model, p, z)
Second virial coefficient | B=second_virial_coefficient(model, T, z)
Mixing | Mixing Function | prop_mix=mixing(model,p,T,z,property)
Excess | Excess Function | prop_E=excess(model,p,T,z,property)
VLE (single-component) | Saturation pressure | p_sat, V_l, V_v=saturation_pressure(model, T)
Saturation Temperature | T_sat, V_l, V_v=saturation_temperature(model, p)
Enthalpy of Vapourisation | H_vap = enthalpy_vap(model, T)
Critical point (single-component) | Critical point | Tc, pc, Vc=crit_pure(model)
VLE and VLLE (multi-component) | Bubble pressure | p_bub, V_l, V_v, y=bubble_pressure(model, T, x)
Bubble temperature | T_bub, V_l, V_v, y=bubble_temperature(model, p, x)
Dew pressure | p_dew, V_l, V_v, x=dew_pressure(model, T, y)
Dew temperature | T_dew, V_l, V_v, x=dew_temperature(model, T, y)
LLE pressure | p_LLE, V_l, V_ll, xx=LLE_pressure(model, T, x)
LLE temperature | T_LLE, V_l, V_ll, xx=LLE_temperature(model, p, x)
Azeotrope pressure | p_az, V_l, V_v, x=azeotrope_pressure(model, T)
Azeotrope temperature | T_az, V_l, V_v, x=azeotrope_temperature(model, p)
VLLE pressure | p_VLLE, V_l, V_ll, V_v, x, xx, y=VLLE_pressure(model, T)
VLLE temperature | T_VLLE, V_l, V_ll, V_v, x, xx, y=VLLE_temperature(model, p)
Critical points (multi-component) | Critical point (VLE) | Tc, pc, Vc=crit_mix(model, x)
UCST | p_UCST, V_UCST, x_UCST=UCST_mix(model, T)
UCEP | T_UCEP, p_UCEP, V_l, V_v, x, y=UCEP_mix(model)

###  4.1 Bulk properties

Bulk properties are properties attributed to a single phase at a given
pressure, temperature and composition (as summarised in table 3). We make use
of the automatic differentiation (AD) provided by the ForwardDiff.jl package91
to provide the derivatives needed to obtain these properties, which ensures
that all properties are compatible with all thermodynamic models. The package
works by providing a Dual number type, that propagates the evaluation of a set
of primal functions ($+$, $-$, sin, exp) with it’s derivative. The derivative
of an arbitrary function composed of differentiable primitives is just the
propagation of those derivatives via the chain rule. This paradigm is known as
differentiable programming.92 Because the differentiation is done at the
primitive level, the error associated with the derivative is similar to the
error of the original function. Furthermore, one differentiation evaluation
results in the original function and its derivative, allowing the use of both
values with no additional evaluation time.

An example of how automatic differentiation is used to obtain the chemical
potential is shown below:

⬇

function VT_chemical_potential(model::EoSModel, V, T, z=SA[1.])

fun(x) = eos(model,V,T,x)

return ForwardDiff.gradient(fun,z)

end

It is important to note that these thermodynamic properties are still defined
in terms of $(V,T,\mathbf{z})$. Thus, if we want to obtain properties in at a
specified $(p,T,\mathbf{z})$, a solver needs to be used to make that
transformation. Clapeyron.jl provides two methods to solve for the volume of a
system:

  * •

Cubic-only method: For a cubic equation of state:

$c_{0}+c_{1}V+c_{2}V^{2}+c_{3}V^{3}=0\,,$ (15)

it is very easy to obtain all roots analytically. When there are three real
roots, one may determine which phase is the most stable by simply checking
which of the three roots corresponds to the smallest Gibbs free energy.

  * •

General method: Most other equations available in Clapeyron.jl cannot be
solved as easily as 15. In the general case, we treat this as a non-linear
equation with multiple roots93:

$p_{0}-p(V,T_{0},\mathbf{z}_{0})=0\,.$ (16)

We make use of an iterative root-finding scheme to solve for the volume. We
use a modified Newton-Raphson step, effectively solving for the logarithm of
the volume, primarily to improve the stability in the liquid phase as the
derivatives can be very large.:

$\displaystyle\ln(V_{i+1})=\ln(V_{i})-\Delta\,,$ (17)
$\displaystyle\Delta=(p_{0}-p(V_{i}))\beta_{T}(V_{i})\,,$ (18)

where $\Delta$ is the step size and $\beta_{T}$ is the isothermal
compressibility. The derivative for the step size of the algorithm above is
obtained from automatic differentiation to improve the numerical stability of
this method.

However, the above only describes how one obtains a single solution, without
checking if it is the most stable solution. This is not a problem if we
already know the phase of our system at the given conditions; this can be
specified as an additional input to all bulk properties listed in table 3:

⬇

julia> model = PCSAFT([”water”]);

\parjulia> u = speed_of_sound(model,1e5,373.15;phase=:liquid)

2142.9834287736817

\parjulia> u = speed_of_sound(model,1e5,373.15;phase=:vapour)

510.1968680367978

In order to ensure that we solve for the specified phase, an initial guess
close to the ‘true’ solution is used:

    * –

Vapour phase: Given that the second virial coefficient can be solved for
directly in Clapeyron.jl, it is possible to obtain a second-order Virial
expansion approximation of all equations of state in Clapeyron.jl:

$Z=1+B(T)\rho\,,$ (19)

where $B(T)$ is the second virial coefficient of the system. It is very easy
to obtain the roots to the above expression; this will provide us with an
initial guess for the volume that is very close to the true solution, as long
as the system is not too close to the critical point ($\sim 0.97T_{c}$).

    * –

Liquid phase: Developing a generalised approach to obtaining good initial
guesses for the liquid volume is very difficult. However, what is generally
true for all equations of state supported in Clapeyron.jl is that they all
have a lower bound for their volume, below which the equation is undefined. A
simple example in SAFT-type equations is that they (should) all be undefined
above a packing fraction of one:

$\eta=1=\frac{\pi}{6}\frac{N}{V}\sum_{i}x_{i}m_{i}\sigma_{i}^{3}\,;$ (20)

it is possible to re-arrange the above expression to obtain the lower-bound
for $V$ ($V_{\text{lb}}$). We can generally assume the liquid phase, due to
its incompressibility, will have a volume to within some factor of the lower
bound (i.e., $V_{\text{liq.}}\approx CV_{\text{lb}}$). Using multiple
dispatch, we can provide tuned $C$ parameters for various equations of state,
to give us more-reliable initial guesses for the liquid phase for various
equations of state. For example:

⬇

function x0_volume_liquid(model::SAFTVRMieModel,T,z)

v_lb = lb_volume(model,z)

return v_lb*1.5

end

In the case where we do not know which phase is more stable, Clapeyron.jl will
solve for the volume twice using both initial guesses. Whichever root returns
the smaller Gibbs free energy is taken to be the true root. To speed this
process up, by default, the two roots are solved for in parallel. It is
possible to switch this feature off using an optional argument (threaded) to
all the bulk properties.

However, this method is not perfect. Naturally, one of its weaknesses is that
we assume we provide good initial guesses for the both phases. Whilst we have
found our guesses to be generally reliable for most conditions, for more-
exotic systems (mixtures with a lot of association or extremely size-
asymmetric mixtures), this may not be the case.

Furthermore, some of these equations of state will observe three seemingly
stable roots (GERG-200827 and PC-SAFT, as pointed out by Privat et al. 93, can
observe an additional root). Given Clapeyron.jl only solves for two roots,
there is no guarantee that the third root will not be the ‘true’ root.

Unfortunately, these problems are difficult to address on a global scale.
However, there is no guarantee either that commercial tools will not
experience similar issues. The difference with an open-source package such as
Clapeyron.jl is that, should these issues arise, they can be dealt with by the
users themselves. For Clapeyron.jl specifically, this can be done in a simple
fashion by modifying built-in functions at the front-end. For example, if the
initial guess for the SAFT-VR Mie equation of state is too small, one can
modify the function:

⬇

import Clapeyron: x0_volume_liquid, SAFTVRMieModel, lb_volume

function x0_volume_liquid(model::SAFTVRMieModel,T,z)

v_lb = lb_volume(model,z)

return v_lb*2.0

end

Julia will automatically over-ride the x0_volume_liquid function, instead
using the one provided by the user. This is an important advantage of open-
source packages; as the package is used more often, more issues are raised and
addressed, quickly. We hope, over time, more users will contribute to the
package, gradually improving the robustness of its solvers.

Once we have solved for the volume, all other bulk properties can be obtained
using ForwardDiff.jl. This has already been illustrated in figures 7(a) and
7(b) where, in the case of the latter, we can see the solver automatically
detects the phase change.

However, another type of bulk property we might be interested in is mixing and
excess functions. These can be of use in continuum models such as the Cahn-
Hilliard equation6. As presented in table 3, the mixing and excess functions
are provided using a single function, rather than for each bulk property. This
allows for users to obtain mixing or excess functions for any bulk property
directly, for example:

⬇

julia> model = PCSAFT([”methanol”,”cyclohexane”]);

\parjulia> mixing(model, 1e5, 313.15, [0.5,0.5], volume)

1.004651584989827e-6

These functions are exemplified in figures 8(a) and 8(b).

(a) Excess volume of ethanol (1) + water (2) at 298.15 K and 1.013 bar.94

(b) Excess enthalpy of acetone (1) + water (2) at 278.15 K and 1.013 bar.95

Figure 8: Mixing functions of various properties available within
Clapeyron.jl.

###  4.2 Pure-component VLE and critical properties

For systems with a single species, in term of vapour–liquid equilibria
properties, there are really only two types of points we are interested in
finding:

  * •

Saturation points: Here, one must find the volumes corresponding the vapour
and liquid phases in equilibrium with one another which satisfy the following
equations:

$\displaystyle p(V_{l},T)$ $\displaystyle=p(V_{v},T)\,,$ (21)
$\displaystyle\mu(V_{l},T)$ $\displaystyle=\mu(V_{v},T)\,,$ (22)

where $\mu$ is the chemical potential of a single component and the subscript
$l$ and $v$ denote the liquid and vapour phases, respectively. For all
equations of state of interest, the above equations represent a non-linear
system of equations which needs to be solved using iterative methods36. As
such, initial guesses must be provided. We have developed two strategies to
obtain these initial guesses:

    * –

van der Waals approximation: Based on the approximations provided by Berberan-
Santos et al. 96, it is possible to obtain the vapour and liquid roots for the
van der Waals equation explicitly. If we map our equations of state to the van
der Waals equation (a mapping that requires the second virial coefficient and
the liquid volume at the specified temperature), it is then possible to obtain
good initial guesses for the liquid and vapour volumes.

    * –

‘Vapour-like’ and ‘liquid-like’ volumes: The van der Waals approximation is
not perfect and, as such, we use the approach developed in the previous
section for the volume solver to obtain initial guesses for the liquid phase
and vapour phase. In the case of the vapour phase, we must supply a value for
the pressure which is usually chosen to be $p=-RT/(4B(T))$. This then
corresponds to a vapour volume of:

$V=-2B(T)\,.$ (23)

This approach is usually used when one is close to the critical point.

The combined approaches above provide a very robust saturation-pressure
solver. However, if our initial guesses fail, we allow users to provide their
own initial guesses using an optional argument
(saturation_pressure(model,T,v0)).

We do note that, another advantage of the ForwardDiff.jl package is our
ability to obtain the Jacobian of equations 21 and 22 meaning that derivative-
based methods for solving linear of systems can be used, allowing for improved
stability of the solvers (provided by the NLSolvers.jl package).

However, the above only allows us to solve for the saturation pressure at a
given temperature. In order to solve for the saturation temperature at a given
pressure ($p_{0}$), an additional layer of iterations needs to be added. This
involves solving the following equation for the temperature:

$p_{0}-p_{\text{sat}}(T)=0\,,$ (24)

where $p_{\text{sat}}(T)$ is the saturation pressure at a given temperature.
To solve the above equation, we use a bracketed method provided by the
Roots.jl package; the solution is bounded between 30% and 100% of the critical
temperature so as to ensure there is always a solution for the saturation
pressure (below 30% of the critical temperature usually corresponds to the
solid phase for most species).

  * •

Critical points: The critical point for a single component is defined as:

$\displaystyle\left(\frac{\partial^{2}p}{\partial V^{2}}\right)_{T}$
$\displaystyle=0\,,$ (25) $\displaystyle\left(\frac{\partial^{3}p}{\partial
V^{3}}\right)_{T}$ $\displaystyle=0\,.$ (26)

For cubic and empirical equations of state, the critical point of a species is
usually an input parameter thus making it unnecessary to solve the above
system of equations. However, for SAFT-type equations, this is not the case
and we must solve the above system of non-linear equations. The third-order
derivative required is seldom obtained analytically (let alone the higher
order derivative required to obtain the Jacobian of the system of equations).
Thus, historically, the critical point is either solved for by iteratively
solving for the saturation point until $V_{l}=V_{g}$ or not supported at all
(such as in gPROMS38). In the case of Clapeyron.jl, these derivatives can all
be obtained by multiple passes of forward differentiation. This leads to a
visually satisfying implementation of the objective function for the critical
point:

⬇

function obj_crit(model::EoSModel, F, T_c, V_c)

_1 = one(T_c+V_c)

∂²A∂V², ∂³A∂V³ = ∂²³f(model, V_c, T_c, SA[_1])

F[1] = ∂²A∂V²

F[2] = ∂³A∂V³

return F

end

As we are solving for the critical volume and temperature in the above
equations, we must provide initial guesses. This is done on a case-by-case
basis where the critical volume is defined with respect to the lower-bound
volume and the critical temperature is defined with respect to the temperature
scaling (for SAFT-type equations, this is the self-interaction energy,
$\epsilon/k_{\text{B}}$). As an example, for the general SAFT equation, the
initial guess for the critical point is $T_{c,0}=2\epsilon/k_{\text{B}}$ and
$V_{c,0}=V_{\text{lb}}/0.3$. Quite surprisingly, the methods used to solve the
above system of equations have been found to be quite robust, regardless of
chain length and presence of association. Nevertheless, we do allow for custom
initial guesses to be used in case the method does fail.

(a) Enthalpy of vapourisation of water.82

(b) Critical points of various alkanes.82

Figure 9: Examples of VLE and critical properties for single-component systems
within Clapeyron.jl.

These two solvers in combination allow us to obtain the VLE envelopes of any
species using any equation of state within Clapeyron.jl; this has already been
exemplified in figures 4(a), 4(b) and 5(a). More examples of this are shown in
figure 9. In figure 9(a), the saturation and critical point solvers are used
to obtain the enthalpy of vapourisation of water from 400 K to the critical
point predicted using various equations of state. In figure 9(b), the critical
points for the first 10 alkanes in the homologous series are obtained using
various SAFT-type equations of state; as is typically expected from these
equations of state, they overestimate the critical point, with the exception
of BACK-SAFT.

###  4.3 Multi-component VLE and critical properties

The last set of properties supported are the multi-component VLE,
liquid–liquid equilibrium (LLE), vapour–liquid–liquid equilibrium (VLLE) and
critical properties. As shown in table 3, a diverse range of these are
supported in Clapeyron.jl. Generally, these can be divided into three types of
functions:

  * •

Phase–equilibria points: These are the bubble, LLE, azeotrope and VLLE
functions where two or three phase are in equilibrium with one another,
satisfying the following equations:

$\displaystyle T^{\alpha}$ $\displaystyle=T^{\beta}$
$\displaystyle(=T^{\gamma})\quad$ (27) $\displaystyle p^{\alpha}$
$\displaystyle=p^{\beta}$ $\displaystyle(=p^{\gamma})\quad$ (28)
$\displaystyle\mu_{i}^{\alpha}$ $\displaystyle=\mu_{i}^{\beta}$
$\displaystyle(=\mu_{i}^{\gamma})\quad$ $\displaystyle\forall i\in
1,\cdots,n_{c}$ (29)

where the superscripts $\alpha$, $\beta$ and $\gamma$ denote different phases.
It is very straightforward to set-up the objective function for such a problem
(regardless of which phases or how many are in equilibrium). Because of
automatic differentiation, it is equally as simple to obtain all the high-
order derivatives for these functions which may be needed when obtaining the
roots of the problem using derivative-based methods.

The only real difficulty for these functions is providing sufficiently
accurate initial guesses. However, if we assume that our saturation and
critical point solvers for single component systems are robust, it is also
possible to obtain reliable initial guesses:

    * –

Bubble and dew points: In this situation, to obtain the initial guess, we make
the assumption that the system can be approximated as an ideal mixture,
obeying Raoult’s law:

$py_{i}=p_{\text{sat},i}(T)x_{i}\,,$ (30)

we can easily obtain $p_{\text{sat},i}(T)$ for all species, allowing us to
obtain an initial guess for $y_{i}$ (for bubble-point calculations) or $x_{i}$
(for dew point-calculations).

With the saturation point of all the pure-species obtained, we also have the
vapour and liquid volumes of each species. We can use these volumes to obtain
initial guesses for the volumes of the mixture assuming ideal mixing:

$V^{\alpha}=\sum_{i}x_{i}^{\alpha}V_{i}^{\alpha}\,.$ (31)

There is one limitation to this method: when one of the components is
supercritical. In this situation, we first obtain the critical point for all
species and, if one component is indeed supercritical, we set:

$\displaystyle V_{l}$ $\displaystyle=V_{c}\,,$ (32) $\displaystyle V_{v}$
$\displaystyle=1.2V_{c}\,,$ (33)

and:

$\tilde{p}_{\text{sat},i}=p(V_{c},T)\,.$ (34)

where $\tilde{p}$ denotes that this is a fictitious saturation pressure as $T$
is greater than $T_{c}$ in this case. All of the above gives a rough estimate
of what the ‘liquid-like’ and ‘vapour-like’ phases in the critical region
would be for a particular species. These initial guesses have proven to be
sufficiently robust for most systems of interest (fringe cases such as highly-
asymmetric mixtures may struggle).

    * –

Azeotrope points: This method is compatible with only binary mixtures. Here,
we re-use the same initial point as bubble and dew points but we assume that
the composition of both components in the liquid and vapour phase is 0.5.

    * –

LLE points: For mixtures excluding highly size-asymmetric ones, the LLE
envelope is generally symmetric. Thus, for a given $\mathbf{x}$, we assume the
composition in the second liquid phase is $1-\mathbf{x}$. The initial guesses
for the liquid volumes are then obtained using equation 31 with the volumes
being obtained from the liquid volume of the saturation point for each
species.

    * –

VLLE points: In the current implementation, without custom initial guesses,
Clapeyron.jl can solve only for binary heterogeneous azeotropes where the
initial guess for the composition of the liquid phases is [0.25,0.75] and
[0.75,0.25], respectively, and the initial guess for the composition of the
vapour phase is [0.5,0.5]. The initial guess for the volumes is obtained using
equation 31 (with the same treatment in the case where a single component is
supercritical).

Note that the above is applicable to only equations of state as, in this case,
equation 12 can be used to solve for bubble and dew points in the case of
activity-coefficient models. LLE and VLLE points cannot yet be solved for in
Clapeyron.jl for activity-coefficient models.

Nevertheless, whilst these initial guesses work well in the general, simple
cases there are a few limitations, some of which have already been mentioned
above. One slightly more-obvious limitation is, for the bubble / dew point, at
a given temperature / pressure and composition, it is possible for two points
to correspond to different dew / bubble points. In order to deal with these
cases where the initial guesses are not sufficiently close to the true or
desired solution, all of the corresponding functions can take custom initial
guesses as an optional argument, similar to the saturation pressure function.
This gives users greater control over the property estimation methods, in
contrast to ‘black-box’ packages such as gPROMS38 or ASPEN37.

Like the saturation points of pure components, these initial guesses serve
only to find the pressure at these conditions. To find the temperature, one
needs to add an additional layer of iterations to solve for the temperature
corresponding to a given pressure (in a similar fashion to equation 24).

(a) $pT$ isopleth for a four-component natural gas system.97

(b) $pT$ isopleth for a 12-component natural gas system.98

Figure 10: Examples of $pT$ isopleths obtained for large systems using various
thermodynamic models within Clapeyron.jl.

In combination, these methods can be used to generate a range of phase
diagrams. Many of these have been shown in previous sections in the form of
simple binary $pxy$ or $Txy$ diagrams (see figures 4(c), 4(d), 5(b) and 8(b)).
However, these methods can also be used for larger systems such as multi-
component, $pT$ isopleths as shown in figures 10(a) (a four-component system)
and 10(b) (a 12-component system) for multiple thermodynamic models (to our
knowledge, this is the first time, in literature, SAFT-$\gamma$ Mie has been
used for such large systems). Whilst these systems are liquid-natural gas
systems with not-too-disimilar interactions between species, it is still
impressive that Clapeyron.jl is able to generate these figures in seconds.

(a) $pxy$ diagram for a methanol (1) + ethane (2) mixture at 298.15 K.99

(b) $Txy$ diagram for a methanol (1) + cyclohexane (2) mixture at 1.01
bar.100, 101

Figure 11: Examples of VLLE phase diagrams for the methanol+ethane mixture
obtained using various thermodynamic models within Clapeyron.jl.

Returning to binary systems, $pxxy$ and $Txxy$ diagrams are shown in figures
11(a) and 11(b), with VLE, LLE and VLLE behaviour. This highlights that
Clapeyron.jl is able model complex phase equilibria for multiple thermodynamic
models.

  * •

Critical points: Similar to pure components, these are rarely supported in
commercial packages. Even in packages in which they are supported, the
critical point for mixtures is rarely solved for directly. This is primarily
because of its complex definition102:

$\displaystyle\mathcal{L}(V,T,z)=\det{\mathcal{H}(A(V,T,z)/(RT))}$
$\displaystyle=0\,,$ (35)
$\displaystyle\det{\begin{bmatrix}\mathcal{H}^{N-1}(A(V,T,z))\\\
(\nabla\mathcal{L}(V,T,z))^{T}\end{bmatrix}}$ $\displaystyle=0\,,$ (36)

where $\mathcal{H}$ denotes the Hessian operator, $\nabla$ denotes the
gradient operator and the superscript $N-1$ denotes that we are only using the
first $N-1$ rows of the Hessian. Without thinking about how to obtain
analytical expressions for this system of equations, as one can clearly see,
the derivatives involved in this system are very unusual and require that one
have access up to the third order derivatives with respect to composition
(without even considering the need for even higher order derivatives in root-
finding algorithms). However, once again, ForwardDiff.jl comes to our aid and
is able to provide us with all the derivatives that we need:

⬇

function Obj_crit_mix(model::EoSModel,F,z,V,T)

f(x) = eos(model,V,T,x)

H(x) = ForwardDiff.hessian(f,x)/R/T

L(x) = det(H(x))

dL(x) = ForwardDiff.gradient(L,x)

M(x) = [H(x)[1:end-1,:];transpose(dL(x))]

F[1] = L(z)

F[2] = det(M(z))

return F

end

This definition for the critical point is true for both the critical point of
VLE envelopes and Upper Critical Solution Temperatures (UCST). The only way to
solve for one of these over the other is by using initial guesses:

    * –

crit_mix: This function will solve for the critical point of VLE envelopes at
a given composition. The initial guesses for the critical temperature and
volume is obtained using a geometric and arithmetic mean (scaled by their
compositions) of the critical temperature and volume for the pure components,
respectively.

    * –

UCST_mix: This function will solve for the UCST at a given temperature (for a
binary mixture) where initial guess for the composition is [0.5,0.5] and the
initial guess for the volume is obtained using the strategy described in
section 4.1 for the liquid volume at this composition.

The objective function for both of these properties is exactly the same and,
like the previous properties, it is possible for users to provide their own
initial guesses. Thus, it is technically possible for users to use the
crit_mix function to obtain the UCST at a given composition (the opposite will
also be true) if they provide sufficiently accurate initial guesses.

One disadvantage of using this method when obtaining critical points for
mixtures is that the computational costs will scale as, at minimum,
$\mathcal{O}(N^{5})$. Whilst this is still reasonable for binary mixtures, it
does become very costly for larger mixtures. In these cases, it might be
computational cheaper to simply solve for the point where the bubble and dew
points are the same.

(a) $pT$ projection of ethane+decane mixture.103

(b) $pT$ projection of ethane+carbon dioxide mixture.104

Figure 12: Examples of $pT$ projections obtained using various thermodynamic
models within Clapeyron.jl.

Nevertheless, the flexibility of this implementation is again demonstrated in
figures 12(a) and 12(b) where $pT$ projections for ethane+decane and
ethane+carbon dioxide are shown. One can see the shape of the critical curve
of the mixture doesn’t affect Clapeyron.jl’s ability to model such systems.

  * •

End points: These are points defined by a critical point in equilibrium with
another phase. As such, we re-use equations 35 and 36, and add the additional
conditions for phase equilibria between two phases. Unfortunately, this is
only currently supported for binary mixtures in Clapeyron.jl and can only find
UCEPs for type II mixtures. For initial guesses, it is assumed that the
critical / liquid phase has a composition of [0.5,0.5] and the vapour phase
has a composition of [0.75,0.25]. The volumes are obtained using the same
method as for the bubble point.

For this property, initial guesses become very important and, unless they are
close to the true solution, the root finder will struggle to converge.
Nevertheless, like all other VLE and critical properties, users are still able
to provide custom initial guesses, allowing them to improve the convergence of
the method and obtaining all other end-points.

Figure 13: $pT$ project of methanol+cyclohexane mixture, showing the
saturation curves and critical points for each species, their critical, UCST,
VLLE and azeotrope curves, along with the UCEP.

All of the VLE and critical properties listed in table 3 are visually
illustrated in figure 13. This is a $pT$ projection of the methanol+benzene
mixture, having the saturation curves for both components (ending at their
critical points) connected by the critical curve. The azeotrope line ranges
from the critical curve down to the UCEP where the UCST and VLLE curves all
intersect. For visual clarity, this figure was made using only a single
equation (PC-SAFT) but can be made using any other thermodynamic model with
the supported parameters.

##  5 Future work

Clapeyron.jl is, by its very nature, a living entity in the sense that it is
(and will remain) in continual development. Despite the many features
currently available, there are still some desirable features currently in
development:

  * •

More thermodynamic models: Despite having an extensive list of thermodynamic
models available in Clapeyron.jl (see table 1), there are still many left to
implement from the literature. However, with its flexible framework, both the
developers and users of Clapeyron.jl will be able to implement new models
quickly (in the time taken to write this article, three new models were
added).

Furthermore, whilst Clapeyron.jl can currently model a range of complex
fluids, it can also be extended to other systems such as electrolytes. This is
already in development on the ElectrolyteModels branch of the main repository
and should be merged into the master branch soon.

  * •

Flash routines: All the methods listed in section 4 assumed prior knowledge of
the phase equilibria. However, users may not always have this knowledge. Flash
routines should solve for the ‘true’ phase equilibrium at specified conditions
(depending on the flash type). As such, they can be considered one of the most
important tools in such a package. Development of these routines in
Clapeyron.jl is already underway with two methods available: a rudementary
Rachford–Rice algorithm105 and a differential evolution algorithm106. It is
our intention to soon implement the more-powerful, but more-complex flash
routines developed by Michelsen and Mollerup36 and the HELD algorithm107.

  * •

Parameter estimation: Whilst users can plug-in existing parameters into
Clapeyron.jl, obtaining new ones is equally as important. How to carry-out
parameter estimation for thermodynamic models is a highly non-trivial
issue108, 68, involving degeneracies and being highly dependent on the dataset
used. Development of a parameter estimation suit has begun on the estimation
branch. Once fully implemented, Clapeyron.jl would be the first open-source
package to provide all features normally provided by commercial software such
as gPROMS38 and ASPEN37 for thermodynamic equations of state.

Whilst the ‘permanent’ development team of Clapeyron.jl is limited to the
three authors of this article, it is our hope, with the gradual increase in
usage of Clapeyron.jl that other individuals in the academic community will
share their expertise and help develop Clapeyron.jl.

##  6 Conclusion

Clapeyron.jl is an open-source package for the implementation, application and
development of thermodynamic models with a strong focus on transparency and
extensibility.

To the best of our knowledge, no other open-source package provides as wide a
selection of thermodynamic models as Clapeyron.jl, which ranges from cubics
equations of state, to SAFT-type equations of state, to activity-coefficient
models, and various empirical Helmholtz models. Our implementation of cubic
equations of state allows for a degree of customisability previously unseen in
thermodynamic packages, giving users a simple way to modify their models
straight from the front end. This level of flexibility is extended to the CPA
equation of state, a unique aspect of our implementation. Whilst not as
modular as cubics, a particularly wide range of SAFT-type equations has also
been provided, all of which can be easily extended to provide the large
variety of variants available in literature.

The property-estimation methods provided in Clapeyron.jl are also extensive,
supporting most properties in which typical users might be interested. In this
aspect particularly, automatic differentiation (through ForwardDiff.jl) has
become a vital tool, not only allowing for easy computation of bulk
properties, but also enabling Clapeyron.jl to support other properties,
particularly critical properties of single and multi-component mixtures, that
even commercial packages are unable to support. Whilst these methods may not
be as robust as their commercial counterparts, users have been provided enough
freedom to have complete control over the solvers and methods to compensate
for this.

Unlike other open-source packages for thermodynamic models, Clapeyron.jl is
readily and easily extensible with its multiple dispatch paradigm, as is the
nature of packages written in Julia. Users can introduce their own parameters
using simple-to-read-and-write CSV files, can implement variants of available
models or entirely new ones and, develop or modify property estimation
methods.

Finally, multiple features such as new thermodynamic models (including
electrolyte models), flash routines and parameter-estimation tools are
currently in development. With the open-source nature of Clapeyron.jl, the
development of new features is not limited to what the ‘permanent’ developers
of the package are focused on. Other members of the community can share their
expertise, develop their own features and integrate them with the package for
the rest of the community to use.

As Clapeyron united the laws of Boyle, Charles, Avogadro and Gay-Lussac in the
ideal gas law, leading to the development of all other thermodynamic models,
we hope that Clapeyron.jl will serve as a unifying framework for all
thermodynamic models and a platform for the development of subsequent models.

##  7 Software Availability

Clapeyron.jl is available through GitHub
(https://github.com/ypaul21/Clapeyron.jl) and can be installed using:

⬇

using Pkg

Pkg.add(”Clapeyron”)

The notebook used to generate all the figures in this article can be found in
the supplementary information and on the GitHub repository
(https://github.com/ypaul21/Clapeyron.jl/examples).

##  8 Supplementary Information

The following documents are provided:

  * •

SI.pdf: This document contains additional information on the structure of the
databases used in Clapeyron.jl.

  * •

user_defined_eos_example.ipynb: This notebook contains an example of how users
can define their own equation of state. The example pertains to the work by
Walker et al. 89. We also provide the csv files needed to run these notebooks
(Ab_initio_SAFT.csv and Ab_initio_SAFT_molarmass.csv).

  * •

We also provide the notebooks needed to generate all the figures presented in
this work:

    * –

cubic_eos.ipynb

    * –

SAFT_eos.ipynb

    * –

activity_models.ipynb

    * –

ideal_eos.ipynb

    * –

mixing_functions.ipynb

    * –

pure_vle_properties.ipynb

    * –

multi-component_vle_vlle_lle_crit.ipynb

All of these are also available on the GitHub repository for Clapeyron.jl.

##  9 Acknowledgements

The authors would like to thank Ramesh Putalapattu, Alex Ames, Sakse Dalum and
Thomas Moore for contributing to the package. Very special thanks go to Dr.
Andrew Haslam who initially introduced PJW and HWY to thermodynamics and has
supported this endeavour since its infancy, providing vital feedback
throughout, including the production of this manuscript.

## References

  * Hendriks et al. 2010 Hendriks, E.; Kontogeorgis, G. M.; Dohrn, R.; de Hemptinne, J.-C.; Economou, I. G.; Žilnik, L. F.; Vesovic, V. Industrial Requirements for Thermodynamics and Transport Properties. _Ind. Eng. Chem. Res._ 2010, _49_ , 11131–11141
  * Müller et al. 2020 Müller, S.; González de Castilla, A.; Taeschler, C.; Klein, A.; Smirnova, I. Calculation of thermodynamic equilibria with the predictive electrolyte model COSMO-RS-ES: Improvements for low permittivity systems. _Fluid Phase Equilib._ 2020, _506_ , 112368
  * Walker et al. 2022 Walker, P. J.; Liang, X.; Kontogeorgis, G. M. Importance of the Relative Static Permittivity in electrolyte SAFT-VR Mie Equations of State. _Fluid Phase Equilib._ 2022, _551_ , 113256
  * Gross and Sadowski 2002 Gross, J.; Sadowski, G. Application of the perturbed-chain SAFT equation of state to associating systems. _Ind. Eng. Chem. Res._ 2002, _41_ , 5510–5515
  * Tihic et al. 2008 Tihic, A.; Kontogeorgis, G. M.; Von Solms, N.; Michelsen, M. L. A predictive group-contribution simplified PC-SAFT equation of state: Application to polymer systems. _Ind. Eng. Chem. Res._ 2008, _47_ , 5092–5101
  * Inguva et al. 2021 Inguva, P. K.; Walker, P. J.; Yew, H. W.; Zhu, K.; Haslam, A. J.; Matar, O. K. Continuum-scale modelling of polymer blends using the Cahn–Hilliard equation: transport and thermodynamics. _Soft Matter_ 2021, _17_ , 5645–5665
  * Walsh et al. 1957 Walsh, J. M.; Rice, M. H.; McQueen, R. G.; Yarger, F. L. Shock-wave compressions of twenty-seven metals. equations of state of metals. _Phys. Rev._ 1957, _108_ , 196–216
  * Hutacharoen et al. 2017 Hutacharoen, P.; Dufal, S.; Papaioannou, V.; Shanker, R. M.; Adjiman, C. S.; Jackson, G.; Galindo, A. Predicting the Solvation of Organic Compounds in Aqueous Environments: From Alkanes and Alcohols to Pharmaceuticals. _Ind. Eng. Chem. Res._ 2017, _56_ , 10856–10876
  * Mahmoudabadi and Pazuki 2021 Mahmoudabadi, S. Z.; Pazuki, G. A predictive PC-SAFT EOS based on COSMO for pharmaceutical compounds. _Sci. Rep._ 2021, _11_ , 6405
  * Cameretti and Sadowski 2008 Cameretti, L. F.; Sadowski, G. Modeling of aqueous amino acid and polypeptide solutions with PC-SAFT. _Chem. Eng. Process._ 2008, _47_ , 1018–1025
  * Clapeyron 1834 Clapeyron, E. Mémoire sur la Puissance Motrice de la Chaleur. _J. Ec. Polytech. (Paris)_ 1834, _XIV_ , 153–90
  * van der Waals 1873 van der Waals, J. D. Over de Continuiteit van den Gas- en Vloeistoftoestand. Ph.D. thesis, University of Leiden, 1873
  * Redlich and Kwong 1949 Redlich, O.; Kwong, J. N. S. On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of Gaseous Solutions. _Chem. Rev._ 1949, _44_ , 233–244
  * Soave 1972 Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. _Chem. Eng. Sci._ 1972, _27_ , 1197–1203
  * Peng and Robinson 1976 Peng, D.-Y.; Robinson, D. B. A New Two-Constant Equation of State. _Ind. Eng. Chem. Fundam._ 1976, _15_ , 59–64
  * Rackett 1970 Rackett, H. G. Equation of state for saturated liquids. _J. Chem. Eng. Data_ 1970, _15_ , 514–517
  * Huron and Vidal 1979 Huron, M.-J.; Vidal, J. New mixing rules in simple equations of state for representing vapour-liquid equilibria of strongly non-ideal mixtures. _Fluid Phase Equilib._ 1979, _3_ , 255–271
  * Boston and Mathias 1980 Boston, P.; Mathias, P. Proceedings of the 2nd International Conference on Phase Equilibria and Fluid Properties in the Chemical Process Industries. West Berlin, 1980; pp 823–849
  * Twu et al. 1980 Twu, C. H.; Lee, L. L.; Starling, K. E. Improved analytical representation of argon thermodynamic behavior. _Fluid Phase Equilib._ 1980, _4_ , 35–44
  * Péneloux et al. 1982 Péneloux, A.; Rauzy, E.; Fréze, R. A consistent correction for Redlich-Kwong-Soave volumes. _Fluid Phase Equilib._ 1982, _8_ , 7–23
  * Michelsen 1990 Michelsen, M. L. A modified Huron-Vidal mixing rule for cubic equations of state. _Fluid Phase Equilib._ 1990, _60_ , 213–219
  * Twu et al. 1992 Twu, C. H.; Coon, J. E.; Cunningham, J. R. A new cubic equation of state. _Fluid Phase Equilib._ 1992, _75_ , 65–79
  * Twu et al. 1995 Twu, C. H.; Coon, J. E.; Cunningham, J. R. A new generalized alpha function for a cubic equation of state Part 1. Peng-Robinson equation. _Fluid Phase Equilib._ 1995, _105_ , 49–59
  * Chapman et al. 1990 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
  * Gross and Sadowski 2001 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
  * Wagner and Kretzschmar 2008 Wagner, W.; Kretzschmar, H.-J. In _International Steam Tables: Properties of Water and Steam Based on the Industrial Formulation of IAPWS-IF97_ ; Edition, S., Ed.; Springer, 2008
  * Kunz and Wagner 2012 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_ , 3032–3091
  * Wilson 1964 Wilson, G. M. Vapor-Liquid Equilibrium. XI. A New Expression for the Excess Free Energy of Mixing. _J. Am. Chem. Soc._ 1964, _86_ , 127–130
  * Renon and Prausnitz 1968 Renon, H.; Prausnitz, J. M. Local compositions in thermodynamic excess functions for liquid mixtures. _AIChE J._ 1968, _14_ , 135–144
  * Fredenslund et al. 1977 Fredenslund, A.; Gmehling, J.; Michelsen, M. L.; Rasmussen, P.; Prausnitz, J. M. Computerized Design of Multicomponent Distillation Columns Using the UNIFAC Group Contribution Method for Calculation of Activity Coefficients. _Ind. Eng. Chem. Process Des. Dev._ 1977, _16_ , 450–462
  * Klamt and Schüürmann 1993 Klamt, A.; Schüürmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. _J. Chem. Soc., Perkin Trans. 2_ 1993, 799–805
  * Klamt 1995 Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. _J. Phys. Chem._ 1995, _99_ , 2224–2235
  * Lin and Sandler 2002 Lin, S. T.; Sandler, S. I. A priori phase equilibrium prediction from a segment contribution solvation model. _Ind. Eng. Chem. Res._ 2002, _41_ , 899–913
  * Holderbaum and Gmehling 1991 Holderbaum, T.; Gmehling, J. PSRK: A Group Contribution Equation of State Based on UNIFAC. _Fluid Phase Equilib._ 1991, _70_ , 251–265
  * Ahlers and Gmehling 2001 Ahlers, J.; Gmehling, J. Development of an universal group contribution equation of state. Prediction of liquid densities for pure compounds with a volume translated Peng-Robinson equation of state. _Fluid Phase Equilib._ 2001, _191_ , 177–188
  * Michelsen and Mollerup 2007 Michelsen, M.; Mollerup, J. _Thermodynamic Modelling: Fundamentals and Computational Aspects_ ; Tie-Line Publications, 2007; p 382
  * 201 2016 Aspen Plus. 2016; https://www.aspentech.com/en/%0Aproducts/engineering/aspen-plus
  * 202 2020 gPROMS. 2020; www.psenterprise.com/products/gproms
  * Bell et al. 2014 Bell, I. H.; Wronski, J.; Quoilin, S.; Lemort, V. Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp. _Ind. Eng. Chem. Res._ 2014, _53_ , 2498–2508
  * 40 Bell, C.; Contributors, Thermo: Chemical properties component of Chemical Engineering Design Library (ChEDL). https://github.com/CalebBell/thermo
  * Chaparro and Mejía 2020 Chaparro, G.; Mejía, A. Phasepy: A Python based framework for fluid phase equilibria and interfacial properties computation. _J. Comput. Chem._ 2020, _41_ , 2504–2526
  * Wilhelmsen et al. 2017 Wilhelmsen, Ø.; Aasen, A.; Skaugen, G.; Aursand, P.; Austegard, A.; Aursand, E.; Gjennestad, M. A.; Lund, H.; Linga, G.; Hammer, M. Thermodynamic Modeling with Equations of State: Present Challenges with Established Methods. _Ind. Eng. Chem. Res._ 2017, _56_ , 3503–3515
  * Rehner and Bauer 2022 Rehner, P.; Bauer, G. FeOs - A Framework for Equations of State and Classical Density Functional Theory. 2022; https://github.com/feos-org/feos
  * Bell et al. 0 Bell, I. H.; Deiters, U. K.; Leal, A. M. M. Implementing an Equation of State Without Derivatives: teqp. _Ind. Eng. Chem. Res._ 0,
  * Rehner and Bauer 2021 Rehner, P.; Bauer, G. Application of Generalized (Hyper-) Dual Numbers in Equation of State Modeling. _Front. Chem. Eng._ 2021, _3_
  * Leal 2018 Leal, A. M. M. autodiff, a modern, fast and expressive {C++} library for automatic differentiation. 2018; https://github.com/autodiff/autodiff/
  * Papaioannou et al. 2014 Papaioannou, V.; Lafitte, T.; Avendaño, C.; Adjiman, C. S.; Jackson, G.; Müller, E. A.; Galindo, A. Group contribution methodology based on the statistical associating fluid theory for heteronuclear molecules formed from Mie segments. _J. Chem. Phys._ 2014, _140_ , 054107
  * Weidlich and Gmehling 1987 Weidlich, U.; Gmehling, J. A modified UNIFAC model. 1. Prediction of VLE, hE, and .gamma..infin. _Ind. Eng. Chem. Res._ 1987, _26_ , 1372–1381
  * Walker and Haslam 2020 Walker, P. J.; Haslam, A. J. A New Predictive Group-Contribution Ideal-Heat-Capacity Model and Its Influence on Second-Derivative Properties Calculated Using a Free-Energy Equation of State. _J. Chem. Eng. Data_ 2020, _65_ , 5809–5829
  * JOBACK and REID 1987 JOBACK, K.; REID, R. ESTIMATION OF PURE-COMPONENT PROPERTIES FROM GROUP-CONTRIBUTIONS. _Chem. Eng. Commun_ 1987, _57_ , 233–243
  * Horstmann et al. 2005 Horstmann, S.; Jabłoniec, A.; Krafczyk, J.; Fischer, K.; Gmehling, J. PSRK group contribution equation of state: Comprehensive revision and extension IV, including critical constants and $\alpha$-function parameters for 1000 components. _Fluid Phase Equilib._ 2005, _227_ , 157–164
  * Robinson and Peng 1978 Robinson, D. B.; Peng, D.-Y. _The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs_ ; Gas Processors Association: Tulsa, 1978
  * Voutsas et al. 2004 Voutsas, E.; Magoulas, K.; Tassios, D. Universal mixing rule for cubic equations of state applicable to symmetric and asymmetric systems: Results with the Peng-Robinson equation of state. _Ind. Eng. Chem. Res._ 2004, _43_ , 6238–6246
  * Chapman et al. 1989 Chapman, W.; Gubbins, K.; Jackson, G.; Radosz, M. SAFT: Equation-of-state solution model for associating fluids. _Fluid Phase Equilib._ 1989, _52_ , 31–38
  * Huang and Radosz 1990 Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules. _Ind. Eng. Chem. Res._ 1990, _29_ , 2284–2294
  * Huang and Radosz 1991 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
  * Fu and Sandler 1995 Fu, Y. H.; Sandler, S. I. A Simplified SAFT Equation of State for Associating Compounds and Mixtures. _Ind. Eng. Chem. Res._ 1995, _34_ , 1897–1909
  * Chen and Mi 2001 Chen, J.; Mi, J. G. Equation of state extended from SAFT with improved results for non-polar fluids across the critical point. _Fluid Phase Equilib._ 2001, _186_ , 165–184
  * Kraska and Gubbins 1996 Kraska, T.; Gubbins, K. E. Phase Equilibria Calculations with a Modified SAFT Equation of State. 1. Pure Alkanes, Alkanols, and Water. _Ind. Eng. Chem. Res._ 1996, _35_ , 4727–4737
  * Kraska and Gubbins 1996 Kraska, T.; Gubbins, K. E. Phase Equilibria Calculations with a Modified SAFT Equation of State. 2. Binary Mixtures of n-Alkanes, 1-Alkanols, and Water. _Ind. Eng. Chem. Res._ 1996, _35_ , 4738–4746
  * Blas and Vega 1997 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
  * Kontogeorgis et al. 1996 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
  * Yakoumis et al. 1997 Yakoumis, I. V.; Kontogeorgis, G. M.; Voutsas, E. C.; Tassios, D. P. Vapor-liquid equilibria for alcohol/hydrocarbon systems using the CPA Equation of State. _Fluid Phase Equilib._ 1997, _130_ , 31–47
  * Kontogeorgis et al. 1999 Kontogeorgis, G. M.; Yakoumis, I. V.; Meijer, H.; Hendriks, E.; Moorwood, T. Multicomponent phase equilibrium calculations for water-methanol-alkane mixtures. _Fluid Phase Equilib._ 1999, _158-160_ , 201–209
  * Gil-Villegas et al. 1997 Gil-Villegas, A.; Galindo, A.; Whitehead, P. J.; Mills, S. J.; Jackson, G.; Burgess, A. N. Statistical associating fluid theory for chain molecules with attractive potentials of variable range. _J. Chem. Phys._ 1997, _106_ , 4168–4186
  * Von Solms et al. 2003 Von Solms, N.; Michelsen, M. L.; Kontogeorgis, G. M. Computational and physical performance of a modified PC-SAFT equation of state for highly asymmetric and associating mixtures. _Ind. Eng. Chem. Res._ 2003, _42_ , 1098–1105
  * Lafitte et al. 2013 Lafitte, T.; Apostolakou, A.; Avendaño, C.; Galindo, A.; Adjiman, C. S.; Müller, E. A.; Jackson, G. Accurate statistical associating fluid theory for chain molecules formed from Mie segments. _J. Chem. Phys._ 2013, _139_
  * Dufal et al. 2015 Dufal, S.; Lafitte, T.; Haslam, A. J.; Galindo, A.; Clark, G. N.; Vega, C.; Jackson, G. The A in SAFT: developing the contribution of association to the Helmholtz free energy within a Wertheim TPT1 treatment of generic Mie fluids. _Mol. Phys._ 2015, _113_ , 948–984
  * Aasen et al. 2019 Aasen, A.; Hammer, M.; Ervik, Å.; Müller, E. A.; Wilhelmsen, Ø. Equation of state and force fields for Feynman-Hibbs-corrected Mie fluids. I. Application to pure helium, neon, hydrogen, and deuterium. _J. Chem. Phys._ 2019, _151_ , 064508
  * Papaioannou et al. 2016 Papaioannou, V.; Calado, F.; Lafitte, T.; Dufal, S.; Sadeqzadeh, M.; Jackson, G.; Adjiman, C. S.; Galindo, A. Application of the SAFT-$\gamma$ Mie group contribution equation of state to fluids of relevance to the oil and gas industry. _Fluid Phase Equilib._ 2016, _416_ , 104–119
  * Abrams and Prausnitz 1975 Abrams, D. S.; Prausnitz, J. M. Statistical thermodynamics of liquid mixtures: A new expression for the excess Gibbs energy of partly or completely miscible systems. _AIChE J._ 1975, _21_ , 116–128
  * Hsieh et al. 2010 Hsieh, C. M.; Sandler, S. I.; Lin, S. T. Improvements of COSMO-SAC for vapor-liquid and liquid-liquid equilibrium predictions. _Fluid Phase Equilib._ 2010, _297_ , 90–97
  * Hsieh et al. 2014 Hsieh, C. M.; Lin, S. T.; Vrabec, J. Considering the dispersive interactions in the COSMO-SAC model for more accurate predictions of fluid phase behavior. _Fluid Phase Equilib._ 2014, _367_ , 109–116
  * Reid et al. 1959 Reid, R. C.; Prausnitz, J. M.; Poling, B. E. _The Properties of Gases and Liquids_ , 4th ed.; McGraw-Hill, 1959
  * Lemmon et al. 2009 Lemmon, E. W.; McLinden, M. O.; Wagner, W. Thermodynamic properties of propane. III. A reference equation of state for temperatures from the melting line to 650 K and pressures up to 1000 MPa. _J. Chem. Eng. Data_ 2009, _54_ , 3141–3180
  * Jørstad 1993 Jørstad, O. Equation of State for Hydrocarbons. Ph.D. thesis, University of Trondheim, 1993
  * Magoulas and Tassios 1990 Magoulas, K.; Tassios, D. Thermophysical properties of n-Alkanes from C1 to C20 and their prediction for higher ones. _Fluid Phase Equilib._ 1990, _56_ , 119–140
  * Waals 1890 Waals, J. D. v. d. Molekulartheorie eines Körpers, der aus zwei verschiedenen Stoffen besteht. _Zeitschrift für Physikalische Chemie_ 1890, _5U_ , 133–173
  * Kontogeorgis and Folas 2009 Kontogeorgis, G. M.; Folas, G. K. _Thermodynamic Models for Industrial Applications: From Classical and Advanced Mixing Rules to Association Theories_ ; 2009; pp 1–692
  * Boukouvalas et al. 1994 Boukouvalas, C.; Spiliotis, N.; Coutsikos, P.; Tzouvaras, N.; Tassios, D. Prediction of vapor-liquid equilibrium with the LCVM model: a linear combination of the Vidal and Michelsen mixing rules coupled with the original UNIF. _Fluid Phase Equilib._ 1994, _92_ , 75–106
  * Wong and Sandler 1992 Wong, D. S. H.; Sandler, S. I. A theoretically correct mixing rule for cubic equations of state. _AIChE J._ 1992, _38_ , 671–680
  * Lemmon et al. 2020 Lemmon, E. W.; McLinden, M. O.; Friend, D. G. _NIST Chemistry WebBook, NIST Standard Reference Database Number 69,_ ; National Institute of Standards and Technology: Gaithersburg MD, 2020; Chapter Thermophys
  * Souza et al. 2018 Souza, L. F.; Al Ghafri, S. Z.; Trusler, J. M. Measurement and modelling of the vapor–liquid equilibrium of (CO2 + CO) at temperatures between (218.15 and 302.93) K at pressures up to 15 MPa. _J. Chem. Thermodyn._ 2018, _126_ , 63–73
  * Butcher and Medani 2007 Butcher, K. L.; Medani, M. S. Thermodynamic properties of methanol-benzene mixtures at elevated temperatures. _J. Appl. Chem._ 2007, _18_ , 100–107
  * Bamberger et al. 2000 Bamberger, A.; Sieder, G.; Maurer, G. High-pressure (vapor+liquid) equilibrium in binary mixtures of (carbon dioxide+water or acetic acid) at temperatures from 313 to 353 K. _J. Supercrit. Fluids_ 2000, _17_ , 97–110
  * Barr-David and Dodge 1959 Barr-David, F.; Dodge, B. F. Vapor-Liquid Equilibrium at High Pressures. The Systems Ethanol-Water and 2-Propanol-Water. _J. Chem. Eng. Data_ 1959, _4_ , 107–121
  * Stabnikov et al. 1972 Stabnikov, V.; Matyushev, B.; Protsyuk, T.; Yushchenko, N. On the Equilibrium in the System Ethyl Alcohol - Water at Atmospheric Pressure. _Pishch. Promst. Kiev_ 1972, _15_ , 49–56
  * Mollerup 1998 Mollerup, J. Unification of the two-parameter equation of state and the principle of corresponding states. _Fluid Phase Equilib._ 1998, _148_ , 1–19
  * 89 Walker, P. J.; Zhao, T.; Haslam, A. J.; Jackson, G. Ab initio development of generalized Lennard-Jones (Mie) force fields for predictions of thermodynamic properties in advanced molecular-based SAFT equations of state. _J. Chem. Phys. (In preparation)_
  * Keller and Contributors 2022 Keller, A.; Contributors, Unitful.jl. 2022; https://github.com/PainterQubits/Unitful.jl
  * Revels et al. 2016 Revels, J.; Lubin, M.; Papamarkou, T. Forward-Mode Automatic Differentiation in Julia. 2016,
  * Ma et al. 2018 Ma, Y.; Dixit, V.; Innes, M.; Guo, X.; Rackauckas, C. A Comparison of Automatic Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation Solutions. 2018,
  * Privat et al. 2010 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
  * Ott et al. 1993 Ott, J.; Sipowska, J.; Gruszkiewicz, M.; Woolley, A. Excess volumes for (ethanol+water) at the temperatures (298.15 and 348.15) K and pressures (0.4, 5, and 15) MPa and at the temperature 323.15 K and pressures (5 and 15) MPa. _J. Chem. Thermodyn._ 1993, _25_ , 307–318
  * French 1989 French, H. Excess enthalpies of (acetone + water) at 278.15, 288.15, 298.15, 308.15, 318.15, and 323.15 K. _J. Chem. Thermodyn._ 1989, _21_ , 801–809
  * Berberan-Santos et al. 2008 Berberan-Santos, M. N.; Bodunov, E. N.; Pogliani, L. The van der Waals equation: Analytical and approximate solutions. _J. Math. Chem._ 2008, _43_ , 1437–1457
  * Avila et al. 2002 Avila, S.; Blanco, S. T.; Velasco, I.; Rauzy, E.; Otín, S. Thermodynamic Properties of Synthetic Natural Gases. 1. Dew-Point Curves of Synthetic Natural Gases and Their Mixtures with Water and Methanol. Measurement and Correlation. _Ind. Eng. Chem. Res._ 2002, _41_ , 3714–3721
  * Mørch et al. 2006 Mørch, Ø.; Nasrifar, K.; Bolland, O.; Solbraa, E.; Fredheim, A.; Gjertsen, L. Measurement and modeling of hydrocarbon dew points for five synthetic natural gas mixtures. _Fluid Phase Equilib._ 2006, _239_ , 138–145
  * Ishihara et al. 1998 Ishihara, K.; Tanaka, H.; Kato, M. Phase equilibrium properties of ethane+methanol system at 298.15 K. _Fluid Phase Equilib._ 1998, _144_ , 131–136
  * Jones and Amstell 1930 Jones, D. C.; Amstell, S. The critical solution temperature of the system methyl alcohol–cyclohexane as a means of detecting and estimating water in methyl alcohol. _J. Chem. Soc._ 1930, 1316–1323
  * Marinichev and Susarev 1965 Marinichev, A.; Susarev, M. No Title. _J. Appl. Chem. USSR_ 1965, _38_ , 1582–1584
  * Bell and Jäger 2017 Bell, I. H.; Jäger, A. Calculation of critical points from Helmholtz-energy-explicit mixture models. _Fluid Phase Equilib._ 2017, _433_ , 159–173
  * Reamer and Sage 1962 Reamer, H. H.; Sage, B. H. Phase Equilibria in Hydrocarbon Systems. Volumetric and Phase Behavior of the Ethane-n-Decane System. _J. Chem. Eng. Data_ 1962, _7_ , 161–168
  * Horstmann et al. 2000 Horstmann, S.; Fischer, K.; Gmehling, J.; Kolář, P. Experimental determination of the critical line for (carbon dioxide + ethane) and calculation of various thermodynamic properties for (carbon dioxide +n-alkane) using the PSRK model. _J. Chem. Thermodyn._ 2000, _32_ , 451–464
  * Rachford and Rice 1952 Rachford, H.; Rice, J. Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. _J. Pet. Technol._ 1952, _4_ , 19–3
  * Bonilla-Petriciolet et al. 2017 Bonilla-Petriciolet, A.; Rangaiah, G. P.; Segovia-Hernández, J. G.; Jaime-Leal, J. E. Phase Stability and Equilibrium Calculations in Reactive Systems using Differential Evolution and Tabu Search. 2017; http://www.worldscientific.com/doi/abs/10.1142/9789813207523_0007
  * Pereira et al. 2012 Pereira, F. E.; Jackson, G.; Galindo, A.; Adjiman, C. S. The HELD algorithm for multicomponent, multiphase equilibrium calculations with generic equations of state. _Comp. Chem. Eng._ 2012, _36_ , 99–118
  * Dufal et al. 2014 Dufal, S.; Papaioannou, V.; Sadeqzadeh, M.; Pogiatzis, T.; Chremos, A.; Adjiman, C. S.; Jackson, G.; Galindo, A. Prediction of Thermodynamic Properties and Phase Behavior of Fluids and Mixtures with the SAFT-$\gamma$ Mie Group-Contribution Equation of State. _J. Chem. Eng. Data_ 2014, _59_ , 3272–3288