ipqa-research/yaeos

Thermodynamic Equations of State, Fortran library with both automatic and anallytical derivation capabilities

There are multiple open source equations of state libraries, like:

• Clapeyron julia
• FeOs rust with Python bindings
• teqp C++ with Python bindings
• thermo python
• thermopack Fortran with Python bindings
• CoolProp C++ with Python bindings

Here we are presenting yet another (still in development) one, that tackles the same problem just, in another way. Mostly exploiting the readability and extensibility of Modern Fortran for scientists to have an easy way to implement new thermodynamic models without dealing with lower-level languages but still getting decent performance. And also this framework provides the possibility of using analytically obtained derivatives so both options are easily available.

This is an experimental work in progress and we recommend the before mentioned libraries if you are intending to use some of this in real work. Big part of the code comes from a refactoring process of older codes so not all parts are easily readable, yet.

We focus mainly on that the addition of a new thermodynamic model as easily as possible. Also providing our models too!

For now, we only include residual Helmholtz model (like Cubic or Saft Equations of State). But we'll be adding other models like $G^E$ (UNIFAC for example).

A little taste of yaeos

A lot of users get the bad picture of Fortran being old and archaic since most of the codes they've seen are written in ancient F77.

use yaeos, only: PengRobinson76, pressure, ArModel

integer, parameter :: n=2   ! Number of components
real(8) :: V, T, P, dPdN(n) ! variables to calculate
class(ArModel), allocatable :: model ! Model

real(pr) :: z(n), tc(n), pc(n), w(n), kij(n, n), lij(n, n)

z = [0.3, 0.7]
tc = [190., 310.]
pc = [14., 30.]
w = [0.001, 0.03]
kij = reshape([0., 0.1, 0.1, 0.], [n,n])
lij = kij / 2 

model =  PengRobinson76(tc, pc, w, kij, lij)

V = 1
T = 150

call pressure(model, z, V, T, P)
print *, P

! Obtain derivatives adding them as optional arguments! 
call pressure(model, z, V, T, P, dPdN=dPdN)
print *, dPdN

How to install/run it

yaeos is intended to use as a fpm

You can either:

• Generate a new project that uses yaeos as a dependency

fpm new my_project

In the fpm.toml file add:

[dependencies]
yaeos = {git="https://github.com/ipqa-research/yaeos"}

• Clone this repository and just modify the executables in the app directory

git clone https://github.com/ipqa-research/yaeos
cd yaeos
fpm run

Including new models with Automatic Differentiation.

We are using the hyperdual module developed by Philipp Rehner and Gernot Bauer

The automatic differentiation API isn't fully optimized yet so performance is much slower than it should be.

type, extends(ArModelAdiff) :: YourNewModel
   type(Substances) :: composition
   real(8), allocatable :: parameters(:) ! A vector of parameters
contains
   procedure :: Ar => arfun
   procedure :: get_v0 => v0
end type

subroutine arfun(self, n, v, t, Ar)
   class(YourNewModel), intent(in) :: self
   type(hyperdual), intent(in) :: n(:) ! Number of moles
   type(hyperdual), intent(in) :: v ! Volume [L]
   type(hyperdual), intent(in) :: t ! Temperature [K]
   type(hyperdual), intent(out) :: ar_value ! Residual Helmholtz Energy

   ! A very complicated residual helmholtz function of a mixture
   Ar = sum(n) * v * t / self%parameters(1)
end subroutine

function v0(self, n, p, t)
   class(YourNewModel), intent(in) :: self
   real(pr), intent(in) :: n(:) ! Number of moles
   real(pr), intent(in) :: p ! Pressure [bar]
   real(pr), intent(in) :: t ! Temperature [K]
   real(pr) :: v0

   v0 = self%parameters(3)
end function

A complete implementation of the PR76 Equation of State can me found in example/adiff/adiff_pr76.f90.

All the thermodynamic properties can be found in yaeos_thermoprops and called like:

use yaeos_thermoprops, only: fugacity_vt
use my_new_model, only: YourNewModel
...
type(YourNewModel) :: eos
eos%parameters = [1, 2, 3]
call fugacity_vt(eos, n, v, t, lnfug=lnfug, dlnphidn=dlnphidn)

Documentation

The latest API documentation for the main branch can be found here. This was generated from the source code using FORD. We're working in extending it more.

yaeos's Issues

Add a consistency check for alfa functions.

This paper seemed interesting, could be a nice addition to generate a tester function for alpha functions

https://www.sciencedirect.com/science/article/pii/S0378381216303648

UNIFAC

Ya hay comenzada una implementación de UNIFAC, pero no está dando buenos resultados y aún faltan derivadas

Federico Benelli: https://github.com/ipqa-research/yaeos/blob/dev/src/models/excess_gibbs/group_contribution/unifac.f90#L138 << Esa línea está complicada por algún motivo
yaeos/src/models/excess_gibbs/group_contribution/unifac.f90 at dev · ipqa-research/yaeos

Federico Benelli: O sea, esa ecuación
Federico Benelli: https://thermo.readthedocs.io/thermo.unifac.html#thermo.unifac.UNIFAC.lngammas_r << Referencia

Phase equilibria procedures

Include phase equilibria related procedures.

• Flash PT/VT
• Phase envelopes
• etc

BUG: Flash doesn't really provide the true light phase

Right now the flash routine assures that the y fraction is that of the lighter phase. While it generally works, the approach given now in the flash routine is not the most robust one, since it compares molar volumes and the true density of the phase can differ a lot (or even invert concerning molar volumes) on highly asymmetric mixtures.

For now I would just remove that option and just leave the user to know which phase is lighter. A feasible way to do this would be to convert to mass volumes but I don't like the idea of forcing the user to provide molar weights just for this.

Add an API to use external libraries

An API to use other libraries would be beneficial to use more models, also for comparisons.

First libraries in mind:

• teqp
• FeOs

Cubic d1 and d2 mixing rules

While there are no much real cases where $\delta_1$ and $\delta_2$ have a mixing rule, including them in the cubic function makes sense

Ideal models

A lot of thermodynamic properties require an ideal model, a new Class to handle this should be provided. And probably there should later be a higher order class that include the different kind of Objects (Ar model, Ideal model, etc)

Add new models

Right now we only provide the PengRobinson76 EoS both with analytical and automatic derivatives. We need to extend this more to more models

Step setter for dS in continuation method

Python API

A Python API for accessibility should be added. Following a singleton pattern in an isolated module could be a good choice

Define some kind of IO to work with besides minimal Fortran scripts

While the library has reached a mature enough API to be straightforward to work with using minimal Fortran scripts (running them with fpm) a good IO with input files can be helpful to provide also sharable executable(s) that do reestablished calculations.

My first thought is to use some standard like json since it makes more sense to use a preestablished format that can also provide some interoperability between programs. Other alternatives are using a scriptable and embeddable language like lua that could give more flexibility.

The ultimate input file will probably be using the Python interface already being worked on its branch.

Work on the documentation

There is little to none documentation since the API rework. So it should be updated.

Also user level pages should be added to ease up the understanding process for the end-user.

BUG: Flash routine gives an error when the iterations variable isn't provided

There is no check of presence on the iterations count variable, this should be an easy fix using an internal counter and set the iterations variable if provided.

Add a consistency module to test any Ar model's Ar and fugacity derivatives

Create a module to run all consistency tests listed in Michelsen and Mollerup chapter 2.3 for any instantiated Ar model.

Gibbs excess models????

Hi

Do you plan to add some Gibbs excess models?

For example, UNIFAC and then build cubic models with Huron-Vidal mixing rules.

I see PSRK, i happy 🐒