Two roots for Carbon Dioxide at standard conditions?

[jonnymaserati]

Hi there

Any idea why the following gives a root in the liquid phase? I'd like to use this for determining the phase (and density) for differing P and T, but seems that it's not able?

import numpy as np
from pint import UnitRegistry
import pyforfluids as pff

model = pff.models.GERG2008()
ureg = UnitRegistry()
Q_ = ureg.Quantity

fluid = pff.Fluid(model, {"carbon_dioxide": 1.0}, temperature=Q_(0, ureg.celsius).to(ureg.kelvin).m, pressure=Q_('1 atm').to('Pa').m)

The last line gives a UserWarning: Two roots were found! Vapor-phase value will be used.

Ref the diagram below, 0 Celsius and 1 atm looks to me like it's well within the gas phase, so why is a fluid phase root being determined?

>>> fluid.density_iterator(Q_('1 atm').to('Pa').m)
(20.16256859062894, 0.044918356427274264, False)

[fedebenelli]

Hello! Thanks for using the package!

This two root problem could be related to the model itself, the GERG equation of state has some weird roots sometimes (still I don't find it logic for a system so simple as this one)

Could you draw an isotherm to better show what the model predicts?

import numpy as np
import matplotlib.pyplot as plt

densities = np.linspace(0.001, 25, 1000)
isotherm = fluid.isotherm(densities)

plt.plot(isotherm["density"], isotherm["pressure"])

[jonnymaserati]

Thanks for making the package and for the quick response!

As requested...

[fedebenelli]

I see, there you can see that there are two roots (in the model world, not necessarily the real world) at those conditions. I'm not sure if the GERG model accurately predicts pure carbon dioxide properties (it's a model developed for natural gas - high in CH4 - mixtures). I'll give this a look later and give you a better hindsight :)

[jonnymaserati]

Looks like there's 5?

[fedebenelli]

Those are unstable roots, part of the mathematical model (like with classic cubic eos where you have one unstable root under critical temperatures)

[jonnymaserati]

Thank you. Is the process then to work inward from the ends and find the first root and stop and if there's a vapor root to use that in preference of a liquid?

[fedebenelli]

The more robust way is to get both roots and calculate either the Gibbs energy or fugacity at those roots (you can get both roots and set the different densities with density_iterator as you've shown. You can check Gibbs energy with fluid["gibbs_energy"]. Or just get the root that you want (either liquid or vapor) and use that

[jonnymaserati]

Hi @fedebenelli

With the following code:

fluid = pff.Fluid(model, {"carbon_dioxide": 1.0}, temperature=Q_(10, ureg.celsius).to(ureg.kelvin).m, pressure=Q_('1 atm').to('Pa').m)
liquid_density, vapor_density, single_phase = fluid.density_iterator(Q_('1 atm').to('Pa').m)

and then:

>>> fluid.set_density(liquid_density)
>>> fluid.calculate_properties()
>>> print(f"liquid phase:\n{fluid.density = }, {fluid.properties.get('gibbs_free_energy') = }")
liquid phase:
fluid.density = nan, fluid.properties.get('gibbs_free_energy') = nan

followed by:

>>> fluid.set_density(vapor_density)
>>> fluid.calculate_properties()
>>> print(f"vapor phase:\n{fluid.density = }, {fluid.properties.get('gibbs_free_energy') = }")
vapor phase:
fluid.density = 0.043298915227438006, fluid.properties.get('gibbs_free_energy') = -5930.668659357086

It seems that the result is not valid if the liquid_density is utilized in the above example - so can't an internal check be done to determine this to prevent an invalid result being returned for the density - I don't think the user should be left to have to check for this?

Is there any way to return the state for a given condition (vapor, liquid, solid or super critical)?

[fedebenelli]

Sorry, right now we don't have an option to get a desired root by default (I agree in that it should be integrated). As a temporary fix, you could use a function like:

def get_liquid(fluid):
    density_liquid = fluid.density_iterator(fluid.pressure)[0]
    fluid_liquid = fluid.copy(density=density_liquid)
    return fluid_liquid

def get_stable(fluid):
    density_liquid, density_vapor, single_phase = fluid.density_iterator(fluid.pressure)
    if not single_phase:
        fluid_liquid = fluid.copy(density=density_liquid)
        fluid_vapor = fluid,copy(density_density_vapor)
        g_v, g_l = fluid_vapor["gibbs_free_energy], fluid_liquid["gibbs_free_energy"]
        if g_v < g_l and not np.isnan(g_v):
            return fluid_vapor
        elif g_l < g_v and not np.isnan(g_l):
             return fluid_liquid
    else:
        return fluid

This whole package is on a whole rework how the models work so we aren't focusing on the user API (for now). But in better ways to include more models. But that can be a little hack for you for now.

Thanks again for using the package and giving some hindsight!

[jonnymaserati]

Thanks @fedebenelli!

Would be nice to include this in the class... I'll see if I can fork it, change it and then PR it?

[fedebenelli]

Sure thing! Any contribution is welcome :) Just make sure to satisfy our CI, we use tox for that: tox