Building spinodal data?

[nick.ma...@gmail.com]

Hi there, thanks for creating CoolProp! This is a most useful library.

I'm interested in calculating the location of the spinodal for single component fluid systems (and ultimately binary systems). With version 6.0.0, I obtained this along an isotherm (for example) by finding:

 

CP.PropsSI('d(P)/d(Dmass)|T','T',T,'Dmass',Dmass_vector,fluid)

 

and finding the first inflection from each side of the saturation curve. Trapezoidal integration of this derivative from the saturation curve provides the pressure at the spinodal at the given temperature.
Repetition with a range of pressures builds a (P,Dmass,T) diagram of the spinodal.
Clearly, this is an inelegant solution, but provided the data I needed.

I upgraded to version 6.1.0 and tried (unsuccessfully) to use the low-level interface get_spinodal_data and build_spinodal methods. Inspection of the source code for AbstractState.h (lines 383-384, http://www.coolprop.org/_static/doxygen/html/_abstract_state_8h_source.html) indicates that it isn't implemented:

383  virtual void calc_build_spinodal(){ throw NotImplementedError("calc_build_spinodal is not implemented for this backend"); };

384  virtual SpinodalData calc_get_spinodal_data(){ throw NotImplementedError("calc_get_spinodal_data is not implemented for this backend"); };

Am I reading this correctly in interpreting that the get_spinodal_data and build_spinodal methods are currently not implemented? Or is it for a particular type of equation of state, etc.


Thanks!

Nick Mason-Smith

[nick.ma...@gmail.com]

Corrigendum: note that

finding the first inflection from each side of the saturation curve.

should read

finding the first read zero crossing from each side of the saturation curve.

For a bit more context: no problems building the phase envelope (with reference to https://github.com/CoolProp/CoolProp/issues/1544) but the code crashes when I try to build the spinodal. I also switched to HelmholtzEOSMixtureBackend as there isn't the same comment about build_spinodal not being implemented (from http://www.coolprop.org/_static/doxygen/html/_helmholtz_e_o_s_mixture_backend_8cpp_source.html#l03722):

Running the code below fails on the call to build_spinodal:

from CoolProp import CoolProp as CP
from matplotlib import pyplot as plt
import numpy as np

# Phase envelope
HEOS = CP.AbstractState('HelmholtzEOSMixtureBackend','R134a')
HEOS.build_phase_envelope('dummy')
p_e = HEOS.get_phase_envelope_data()
plt.plot(p_e.T,np.divide(p_e.p,1e6))
plt.xlabel('$T$ (K)')
plt.ylabel('$P$ (MPa)')

# Spinodal
HEOS.build_spinodal()

 

Thanks in advance!


Cheers,

[nick.ma...@gmail.com]

Solved my own problem - turns out HelmholtzEOSMixtureBackend actually wants a mixture! Changed to use a mixture (with a trivial amount of the second component) and it worked.

Is the build_spinodal object likely to be added for pure fluids?

A snippet regarding my method - in case anyone is interested in generating plots of the binodal and spinodal in a Dmass-molefrac-T space:

from CoolProp import CoolProp as CP
from matplotlib import pyplot as plt

from mpl_toolkits.mplot3d import Axes3D
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
fluidA = 'n-Decane'
fluidB = 'n-Dodecane'

# Pure fluid A phase envelope
HEOS = CP.AbstractState('HEOS','{}'.format(fluidA))
HEOS.build_phase_envelope('dummy')
binodal = HEOS.get_phase_envelope_data()
molar_mass = HEOS.molar_mass()
ax.plot(np.multiply(binodal.rhomolar_liq,molar_mass),np.zeros(len(binodal.T)),binodal.T,'k')

# Mixtures
HEOS = CP.AbstractState('HelmholtzEOSMixtureBackend','{}&{}'.format(fluidA,fluidB))
tol = 1e-3
fluidB_mole_fractions = np.linspace(tol,1-tol,25)
# Mixture phase envelope and spinodal
for i in range(0,len(fluidB_mole_fractions)):
    fluidB_mole_fraction = fluidB_mole_fractions[i]
    fluidA_mole_fraction = 1-fluidB_mole_fraction
    HEOS.set_mole_fractions([fluidA_mole_fraction,fluidB_mole_fraction])
    molar_mass = HEOS.molar_mass()
    HEOS.build_phase_envelope('dummy')
    binodal = HEOS.get_phase_envelope_data()
    ax.plot(np.multiply(binodal.rhomolar_liq,molar_mass),fluidB_mole_fraction*np.ones(len(binodal.T)),binodal.T,'k',label='Binodal')
    ax.plot(np.multiply(binodal.rhomolar_vap,molar_mass),fluidB_mole_fraction*np.ones(len(binodal.T)),binodal.T,'k')       
    HEOS.build_spinodal()
    spinodal = HEOS.get_spinodal_data()
    ax.plot(np.multiply(spinodal.delta,HEOS.rhomolar_reducing())*molar_mass,fluidB_mole_fraction*np.ones(len(spinodal.tau)),np.multiply(np.power(spinodal.tau,-1),HEOS.T_reducing()),'k--',label='Spinodal')
Axes3D.set_xlabel(ax,xlabel='Density (kg/m$^3$)')
Axes3D.set_ylabel(ax,ylabel='{} mole fraction'.format(fluidB))
Axes3D.set_zlabel(ax,zlabel='$T$ (K)')
handles,labels = ax.get_legend_handles_labels()
ax.legend(handles[0:2],labels[0:2])

Thanks,

[Ian Bell]

Hi Nick,

Very cool stuff! Sorry I wasn't able to reply earlier, but you are right, the spinodal tracer is part of the critical locus tracer from my paper on critical points.  So therefore, only for mixtures, not pure fluids.  I'm not 100% sure it works with multicomponent mixtures either :/ 

You would be absolutely welcome to contribute to improving the spinodal tracer for pure fluids.  My hands are currently tied so I can't do it myself, but I can review code contributions.

Best,

Ian

--
You received this message because you are subscribed to the Google Groups "coolprop-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to coolprop-users+unsubscribe@googlegroups.com.
To post to this group, send email to coolprop-users@googlegroups.com.
Visit this group at https://groups.google.com/group/coolprop-users.
For more options, visit https://groups.google.com/d/optout.

[Matthis Thorade]

Hi,
for pure fluids you might also want to take a look at this notebook (and the paper cited in there):
https://github.com/CoolProp/CoolProp/blob/master/doc/notebooks/Maxwell_Loop.ipynb
https://nbviewer.jupyter.org/github/CoolProp/CoolProp/blob/master/doc/notebooks/Maxwell_Loop.ipynb

https://doi.org/10.1134/S0036024406040030

[Nick Mason-Smith]

Thanks Ian and Matthis! I'll have a look at these in more detail soon. Ian thanks for pointing me towards the paper on critical points. I should have time around the end of March, into April, to make some contributions.

Thanks,

____________________________________________________
Nicholas Mason-Smith

PhD Candidate

Laboratory for Turbulence Research in Aerospace and Combustion

Department of Mechanical and Aerospace Engineering

Monash University Clayton Campus
Clayton VIC 3800
Australia

M: +61 (0)439 595 180
E: nick.ma...@gmail.com
https://scholar.google.com.au/citations?user=PBbVMTEAAAAJ&hl=en
https://www.researchgate.net/profile/Nicholas_Mason-Smith

--
You received this message because you are subscribed to a topic in the Google Groups "coolprop-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/coolprop-users/prZVxyKzH0o/unsubscribe.
To unsubscribe from this group and all its topics, send an email to coolprop-users+unsubscribe@googlegroups.com.
To post to this group, send email to coolpro...@googlegroups.com.

[Julian Restrepo]

Dear all, thanks a lot for Coolprop, it is an amazing tool!.

 

I am trying to calculate the Spinodal line for low pressures < 100 kg/m3 in the vapor region for a methane carbon dioxide mixture, nevertheless I can´t use the command HEOS.build_spinodal() because it can't get the data for low pressures. Attached you will find the curves of the Spinodal data obtained using the code of  Nick Mason-Smith (Thanks for the code), as you can see in the plot, it is not possible to obtain a low pressure Spinodal curve.

 

I have tried to do something with python and the scipy.optimize module to find the place where the partial derivative of the pressure and volume becomes zero, I have obtained some results but I am not sure if the data obtained corresponds to the vapor Spinodal curve, because I don´t know if there are several roots in the region above the binodal curve. Attached will find a plot with the Binodal and Spinodal data (spinodal_binodal.png), obtained with this code, for a mixture of 50% Methane, 50% carbon dioxide.

 

import CoolProp as cp

import scipy.optimize as optimize

import numpy as np

import math

import matplotlib.pyplot as plt

 

 

Tmin =200 #Minimum temperature [K]

x0 = 0.5 #Carbon dioxide molar fraction

F= 'Methane&CarbonDioxide' #Fluid String

F2='HEOS::Methane['+str(1-x0)+']&CarbonDioxide['+str(x0)+']' #Fluid String High Level Interface

"Binodal Line construction using CoolProp"

 

HEOS = cp.AbstractState('HelmholtzEOSMixtureBackend',F)

HEOS.set_mole_fractions([1-x0,x0])

molar_mass = HEOS.molar_mass()

HEOS.build_phase_envelope('dummy')

binodal = HEOS.get_phase_envelope_data()

Tmax=math.floor(max(binodal.T)) # Max temp

n= Tmax-Tmin #Number of data

P = np.zeros((n,1))

psat = np.zeros((n,1))

rho = np.zeros((n,1))

T=np.zeros((n,1))

for dt in range(0,n):

    try:

        psat[dt,0]=cp.CoolProp.PropsSI('P','T',Tmin+dt,'Q',1,F2)*0.6

        T[dt,0]=Tmin+dt

        def spinodal(R):

            P[dt,0]=R

            HEOS=cp.AbstractState("HEOS", F)

            HEOS.set_mole_fractions([1-x0,x0])

            HEOS.update(cp.PT_INPUTS,P[dt,0],T[dt,0])

            rho[dt,0]=HEOS.rhomass()

            return (HEOS.first_partial_deriv(cp.iP, cp.iDmass, cp.iT)/rho[dt,0]**2)

        P[dt,0]=(optimize.fsolve(spinodal,psat[dt,0],xtol=1e-04,maxfev=1000))

    except ValueError as VE:

        print(VE)

 

My questions are:

 

I am doing something wrong with CoolProp? I should obtain the Spinodal curve for low pressures with the HEOS.build_spinodal() command?

 

The roots that I find with fsolve, corresponds to the Spinodal curve?

Thanks a lot for your help!!

Kind regards from Brazil.

Julián.

[Matthis Thorade]

Julian, probably Nick is the perfect person to discuss that!

Once you two have some good code that works, e.g. as Jupyter notebook, it would be great to share that here or as a pull request to the Coolprop repo!

[Ian Bell]

Julian,

You might want to consider reading our paper about finding critical points(https://ws680.nist.gov/publication/get_pdf.cfm?pub_id=920895), and also, look into the code that is used to trace spinodals in CoolProp.  Basically we lock onto the spinodal and trace it around like a bucking bronco, marching along the spinodal and capturing all its crazy behaviors.  It's too bad that trace_spinodal doesn't work in your case as that is by far the simplest solution when it works.  What error did you get with that in particular?

Another idea is to calculate the surface of the minimum eigenvalue of the Hessian matrix of the Helmholtz energy density (where the independent variables are the molar concentrations), and see where that is equal to zero (by a contouring algorithm, for instance), as that is an alternative (and easy to evaluate) solution for the spinodal.  Are you a Python person? If so, I can give you a ready-made (though a bit slow) solution for this problem. It relies on the isochoric tracing work that we did recently (https://onlinelibrary.wiley.com/doi/abs/10.1002/aic.16074).  You'll just need a compiler that can handle C++11, which is nowadays easy to come by.

All the best,

Ian

--

You received this message because you are subscribed to the Google Groups "coolprop-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to coolprop-users+unsubscribe@googlegroups.com.
To post to this group, send email to coolprop-users@googlegroups.com.

[Julián Restrepo]

Dear Ian.

Thanks a lot for your kind answer.

I took a long time to answer because I am studying all the concepts involved in the metastable states for real gas mixtures, thanks for all the info.

In the end, I am trying to find the Wilson line for supersaturated vapour (Homogeneous nucleation for a real gas mixture of carbon dioxide + Methane or Nitrogen), for this reason, is essential to me find the spinodal  line for low densities and later the Maxwell isothermal loop between the spinodal and the saturation state.

I have seen that some people have worked on this topic as Mathias and Nick, which is a good starting point, I will work a little more on this subject and publish my code on a Jupyter notebook as Mathias suggest.

So in this context, I will answer your questions:

1. I didn't find any error from python or CoolProp when I run the command (HEOS.get_spinodal_data()), is just that this command doesn't calculate the Spinodal for low pressures.

2. I am not a Python person, but I am learning :).

So if you could give me any clue about how to solve this problem or share the code about the Hessian matrix of the Helmholtz energy density, I will be very grateful.

Kind regards

Julián

To unsubscribe from this group and stop receiving emails from it, send an email to coolprop-user...@googlegroups.com.
To post to this group, send email to coolpro...@googlegroups.com.

[Ian Bell]

Sorry for the slow reply.  I gave the link to my paper at the end of my last response, which includes the description of how to calculate this, and also includes the code in C++.  It uses coolprop as the backend.

[Julian Camilo Restrepo Lozano]

Thanks a lot Ian! I will start to work with it.

 

Have a good day!