Errors with SRK and PR cubic EOS for mixtures at high temperatures and pressures

[Cory Kinney]

I'm trying  to calculate the temperature and pressure after a reflected shock in a shock tube for a real gas using CoolProp to solve for the following thermodynamic variables used in the algorithm:

• hmass
• rhomass
• cpmass
• isothermal_compressibility
• isobaric_expansion_coefficient

The conditions we expect are often between 1000K and 2000K, with pressures of 100 bar to 300 bar, for various mixtures of hydrocarbons (e.g. methane or natural gas), oxygen, argon, and/or carbon dioxide. 

I have been unable to reliably use CoolProp with either the SRK or PR cubic EOS at even a pressure of 100 bar. I encounter quite a few different errors - most commonly, there are three roots for the cubic; however, there are several other errors I've been getting that I can't make sense of. 

I've included details about the environment I'm running my code in, as well as sample code and the output/errors I've gotten. Hopefully I'm not misunderstanding how to use CoolProp, or what it should be able to do, so as to waste anyone's time, but from my limited experience it seems very unpredictable and therefore unusable for me. If anyone is able to offer any advice or guidance, that would be much appreciated. I can always open an issue on the GitHub, but only if these are actual issues with CoolProp and not an issue with my implementation/understanding.

Thank you for your time.

Environment:

Python 3.8 

CoolProp 6.4.1 

Code Outline (".../test.py"):

import CoolProp as CP

backend = "PR"

T = 1500 # [K]
P = 100e5 # [Pa]
mixture = {

"AR": 1
}

EOS = CP.AbstractState(backend, "&".join(mixture.keys()))
EOS.set_mole_fractions(mixture.values())

EOS.update(CP.PT_INPUTS, P, T)
print(f"Density = {EOS.rhomass():.3f} kg/m^3")

Samples Cases/Outputs/Comments:

For all these cases, only the mixture dictionary and backend string are chagned from the above code outline

AR: 1 (PR)

Traceback (most recent call last):

  File ".../test.py", line 14, in <module>

    EOS.update(CP.PT_INPUTS, P, T)

  File "CoolProp\AbstractState.pyx", line 102, in CoolProp.CoolProp.AbstractState.update

  File "CoolProp\AbstractState.pyx", line 104, in CoolProp.CoolProp.AbstractState.update

ValueError: Cubic has three roots, but phase not imposed and guess density not provided

AR: 1 (SRK)

[same as above]

CO2: 1 (PR)

[same as above]

CO2: 1 (SRK)

[same as above]

Seems like neither Argon or CO2 pure fluids are working at T = 1500 K, P = 100 bar

AR: 1, CO2: 0 (PR)

Traceback (most recent call last):

  File ".../test.py", line 15, in <module>

    EOS.update(CP.PT_INPUTS, P, T)

  File "CoolProp\AbstractState.pyx", line 102, in CoolProp.CoolProp.AbstractState.update

  File "CoolProp\AbstractState.pyx", line 104, in CoolProp.CoolProp.AbstractState.update

ValueError: critical point finding routine found 0 critical points

Somehow adding Carbon Dioxide to the AbstractState but setting its mole fraction to 0 resulted in a different error than when the AbstractState was initialized with only Argon despite not changing the actual composition

AR: 0.5, CO2: 0.5 (PR)

Density = 33.054 kg/m^3

Somehow the 50-50 mixture of Argon and CO2 works for PR

AR: 0.5, CO2: 0.5 (SRK)

Traceback (most recent call last):

  File ".../test.py", line 15, in <module>

    EOS.update(CP.PT_INPUTS, P, T)

  File "CoolProp\AbstractState.pyx", line 102, in CoolProp.CoolProp.AbstractState.update

  File "CoolProp\AbstractState.pyx", line 104, in CoolProp.CoolProp.AbstractState.update

ValueError: Unable to find gaseous density for T: 1500 K, p: 1e+07 Pa

The 50-50 mixture doesn't work with SRK, though

AR: 0.3, CO2: 0.7 (PR)

Traceback (most recent call last):

  File ".../test.py", line 15, in <module>

    EOS.update(CP.PT_INPUTS, P, T)

  File "CoolProp\AbstractState.pyx", line 102, in CoolProp.CoolProp.AbstractState.update

  File "CoolProp\AbstractState.pyx", line 104, in CoolProp.CoolProp.AbstractState.update

ValueError: p is not a valid number

This error doesn't make any sense to me given that the 50-50 mixture with the same T, P, and PR didn't give an error

[Ian Bell]

Have your tried the development version or specifying the phase? 

[Cory Kinney]

I just installed the CoolProp 6.4.2.dev0 package for Python 3.8 x64 from the SourceForge and confirmed that for each of the test cases I tried, the results are the same. I was hesitant to specify the phase because this code would ideally be run automatically after an experiment is complete, with the only inputs being the initial temperature and pressure, mixture composition, and shock velocity. We don't always know the exact conditions to expect before an experiment, and I don't want to hard code in a specific phase when it may or may not be accurate.

[Ian Bell]

What does the p(rho) curve look like for your specified T for each model? 

[Cory Kinney]

I'm just starting to learn what CoolProp can do, so I'm not sure if there is a P(rho) curve graph built-in, but here's what I generated using code as similar as possible to how I was generating the errors, only specifying density as opposed to temperature. This graph is for argon at 1500 K from 10 to 40 kg/m^3, which covers about 30 to 130 bar; everything looks about as I would expect. 

Graphing Code:

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

backends = ["PR", "SRK"]

D = np.linspace(10, 40, 100) # [kg/m^3]

T = 1500 # [K]

mixture = {
    "AR": 1
}

plt.figure()
for backend in backends:

    EOS = CP.AbstractState(backend, "&".join(mixture.keys()))
    EOS.set_mole_fractions(mixture.values())

    P = []
    for Dmass in D:
        EOS.update(CP.DmassT_INPUTS, Dmass, T)
        P.append(EOS.p())

    plt.plot(D, P, label=backend)

plt.ylabel("Pressure [Pa]")
plt.xlabel("Density [kg/m^3]")
plt.legend()
plt.show()

[Cory Kinney]

Regarding specifying the phase, I'm not entirely understanding how that would help in picking the right root. From this graph on the High-Level Interface page on the CoolProp documentation, my understanding is that the phase could be inferred by CoolProp if the critical pressure and temperature are known. I would think that for a cubic EOS, whose parameters depend on critical temperature and pressure, that CoolProp would be able to infer the phase and pick the right root. Is it possible that there multiple roots in a single phase, or is it not actually as straightforward to infer the correct phase? I'm just trying to avoid explicitly specifying a phase when I don't know what the conditions are going to be since the algorithm I have to implement is iterative.

[Ian Bell]

The phase determination is only very simple for pure fluids. For mixtures,  the critical point could in principle be anywhere, or there could be multiple critical points. Even with simple EOS like cubic EOS.  For that reason you may need to specify the phase,  or do your own density rootfinding.

--
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-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/coolprop-users/2c0da750-eba5-472c-9d6f-2e2a018a696dn%40googlegroups.com.

[Cory Kinney]

Do you think in the case that there are multiple roots, using the ideal gas EOS as a guess for the starting point might improve the root finding?

[Ian Bell]

Only if you know a priori that your inputs are gas,  in which case you should specify the phase.

The phase spec only deals with the case where there are multiple solutions, so you might want to specify liquid and specify gas and keep the density root that makes sense for your problem. 

It is possible that there are bugs, I'll try to find time to look into it,  but not for at least a week.   In the meantime,  if you could plot p(rho) for each of your compositions to the covolume density,  that would be informative. 

To view this discussion on the web visit https://groups.google.com/d/msgid/coolprop-users/a9b193cc-b770-459c-9fca-6195863c9223n%40googlegroups.com.