longemen3000 / lavoisiercore.jl

Abstract Definitions and Functions to establish a thermodynamic library in julia

LavoisierCore.jl

Abstract Definitions and Functions to establish a thermodynamic library in julia

What?

LavoisierCore intends to provide a thermodynamic property calculator, using helmholtz models and automatic differenciation to calculate all relevant properties with just one function definition. for example, lets say you are a chemist with a brand-new Helmholtz model, but you don't want to implement all the properties that would transform your new EoS in something useful.

how?

for example, let's say you want to implement a new helmholtz model:

struct MyHelmholtzModel <: AbstractThermoModel
...#here you store your constants
...
end

then you

the package is the definition of the function core_helmholtz(model,v,T,x), with v being the molar volume (m3/mol), T being the temperature (Kelvin), and x being the vector of molar fractions. the rest is obtained via autommatic differenciation.

core_helmholtz(model::MyHelmholtzModel,v,T,x)
...
end

And that's it!, with that, all relevant thermodynamic property's functions will be created, using the power of Julia's multiple dispatch and the powerful tools of automatic differenciation available (ForwardDiff at the moment, a reverse AD tool in the future).

at this moment, this package has implementations of

• IAPWS 95 (formulation of water)
• GERG 2008 (21 natural gas compounds)

and the following properties:

• compressibility factor
• pressure
• internal energy
• enthalpy
• entropy
• isochoric heat capacity
• isobaric heat capacity
• molar volume (from pressure and temperature)
• mass volume (from pressure and temperature)
• molar density (from pressure and temperature)
• mass density (from pressure and temperature)

if you defined your helmholtz equation, then you have access to all those functions.

Usage example:

m = IAPWS95() #this model contains everything, so it doesn't need any variables to be created

# with unitful properties, it handles automatically the conversion beetween molar and mass density, 
#molar and mass especific volume, and the usual temperatures

julia> pressure(m,(322.0)u"kg/m^3",647.0u"K",[1.0])
220.64298608634772 bar

julia> pressure(m,55u"cm^3/mol",(647.0)u"K",[1.0])
220.59360203472647 bar

julia> pressure(m,55u"cm^3/mol",(373.85)u"°C",[1.0])
220.59360203472647 bar

#Big Floats Support
julia> pressure(m,big(322.0)u"kg/m^3",big(647.0)u"K",[1.0])
220.6429860863478813957299832315505132338714379480500058877684060654357069055415 bar

This package is in heavy development, don't dare to even try to use this in production, but feel free to submit your own Helmholtz equations,the long term idea is to implement the DiferencialEquations.jl of the thermodynamic equations

Why?

Aspen is expensive, COOLPROP uses C++, and the thermodynamic papers are seldom implemented for free (for example, the SAFT equations), and because is fun!, why not?

#News ##September 19, 2019 Added a equilibrium solver in Volume-Temperature-mol, based on a unified representation of helmholtz equilibria (). needs testing with consistent EOS

##November 3, 2019 Added a PT two phase solver based on the TREND solver (EXPERIMENTAL) (Gernert el. al, 2014). added an interfase to call a solver:

water_and_gases = GERG2008(:H2O,:N2,:O2)
method = Gernert()
P0 = 1.1u"atm" #it accepts a number (Pa), but in this case, unitful quantities are better
T = 25u"°C"
method = Gernert()
x0 = [0.6,0.2,0.2]
phases = pt_flash(water_and_gases,P0,T0,x0,method)
entropy.(model,phases) #vector of entropies with units
core_entropy.(model,phases) #vector of entropies without units (SI)

##November 24, 2019 I need a lot of corrections on the PT solver, im working on it!. meanwhile, there are two new functions: core_pure_p_flash and core_pure_t_flash, that calculates saturation's temperature and pressure of one component: for example:

water = IAPWS95
phases = core_pure_t_flash(water,275.0,[1.0]) #2°C, near the triple point, in IAPWS results
julia> core_pressure.(water,phases)
2-element Array{Float64,1}:
 698.4511870563939 #Pa
 698.451166699837 #Pa
phases2 = core_pure_p_flash(water,698.451166699837,[1.0])
julia> core_temperature.(phases2)
2-element Array{Float64,1}:
 274.99999999999994
 274.99999999999994

a future idea is to group all those flash function under one method: flash, and dispatch appropriately. Any suggestions are appreciated and encouraged! (in other words,please leave suggestions)

lavoisiercore.jl's Issues

Documentation

TODO:

• Write about defining custom gradients
• Write about function interactions and necessary definitions
• Write About the mol_fraction, mol_number, mass_fraction, mass_number functions
• Write about the HelmholtzPhase struct

Add dependencies to project

Gotta add Sobol, Mapped Arrays, NLopt, Optim, BenchmarkTools to the project otherwise it won't build.

If you set up CI/CD it will help with these kinds of 'bugs'.

in need of a consistent EOS to test phase equilibria solvers

A lot of recent methods that solves phase equilibria take the assumption of a consistent EOS: one where there are only two Maxwell loops on the P-V curve, IAPWS95 and GERG2008. are empirical, non-consistent methods, as they have a third maxwell loop where a false pseudostable phase (with a lower gibbs energy than the true equilibrium). global minimizations (like the HELD method) or newton with linesearch minimizations of the hessian fail to account for that, so a method with expensive volume search evaluations is needed for non-consistent EOS.

Full support for Measurements.jl

I recently discover the power for measurements.jl, but some equations don't work. is not a priority, but it will help a lot to the scientific use of the ecuations with real measurements to have EoS with uncertainty support

apropiate property naming

CamelCase? camelCase (CAPE-OPEN convention) ? with_lines_i_between (julia style)

Add configuration interface

is necessary to add a configuration interface, to define tolerance, max iterations, and other things

Some thoughts (combine with Thermodynamics.jl?)

I'm finished with my heavy work this semester and ready to put some work into thermodynamics in Julia!

First, I really like your core code (the core.jl file, which seems to be the abstract Helmholtz models and the base of the HELD algorithm--I couldn't really tell if the HELD algorithm is finished yet but I know it's quite a bit of work so perhaps not?).

Second, I think it would be good to plan how to organize the repository. I really like the DifferentialEquations.jl model, where they have a bunch of independent packages for different things and the DifferentialEquations.jl module really just loads all those other packages and reexports them. I also like the idea of putting the core phase equilibria algorithm and other abstract things into a Base package (like DiffEqBase.jl). I actually created a ThermodynamicsBase.jl but there's not a lot there. On a similar note, I'm not sure LavoisierCore is the best name for the thermodynamics library. I think there's value in having a thermodynamics library outside of its use for a particular process simulator, and that it name should state that it's a thermodynamics library. In that vein I suggest that LaviosierCore and Thermodynamics.jl be combined into a single repository named Thermodynamics.jl. You can be the owner, especially since I plan to work more on molten salt and actinide chemistry/thermodynamics rather than general things, since that's what I need for my research, but I think it's a better name.

Handling of physical units

I see in the core code that there's some use of units from Unitful.jl. While I like the idea of integration with Unitful.jl, I think I would prefer it to be a non-default mode. It could be something like this:

pvtxstate(2.0, missing, 500.0, [0.3,0.4,0.3]) #defaults to bar, m^3/mol, K, but could be easily changed by a function like set_base_units("psi", "ft^3/mol", "R")

Or you could use

physical_units() pvtxstate(2.0u"bar", missing, 500.0u"K", [0.3,0.4, 0.3])

Actually, the physical_units function may be unnecessary--you could just dispatch on whether you receive physical quantities or ordinary numbers as input.