Can I use feos to plot pressure-density diagram of binary mixture?

[RayleighSun]

Hi, I tried to plot pressure-density diagram of CO2-propane binary mixture with 0.5/0.5 mole fraction, under the help of examples/pcsaft_phase_diagram.ipynb (https://github.com/feos-org/feos/blob/main/examples/pcsaft_phase_diagram.ipynb) but I cannot set mole fraction of mixture and my code is not working. Is there a similar demo code to do so? Thanks!

[g-bauer]

Hey,

there is no utility function like PhaseDiagram for constant composition (edit: at given pressures). You can create phase equilibria for a given vapor or liquid composition using PhaseEquilibrium.dew_point or PhaseEquilibrium.bubble_point, respectively.

Maybe you can use this as starting point?

from feos.si import *
from feos.eos import *
from feos.pcsaft import *
import numpy as np
import matplotlib.pyplot as plt

parameters = PcSaftParameters.from_json(
    ["carbon dioxide", "propane"], 
    "../parameters/pcsaft/esper2023.json"
)
pcsaft = EquationOfState.pcsaft(parameters)

x = np.array([0.5, 0.5])

# pressures between 1 BAR and 10 BAR
pressures = SIArray1.linspace(BAR, 10*BAR, 11)
states = []

for pressure in pressures:
    states.append(PhaseEquilibrium.bubble_point(
        pcsaft, 
        pressure, 
        liquid_molefracs=x,  
        tp_init=300*KELVIN)
    )

densities = [
    # you can access the liquid and vapor states of a
    # PhaseEquilibrium object via .liquid and .vapor
    state.liquid.mass_density() / (KILOGRAM / METER**3) 
    for state in states
]

# create plot
plt.plot(pressures / BAR, densities)

[RayleighSun]

Yes, I tried to draw a saturation curve of this binary mixture, is my code right?

from feos.si import *
from feos.pcsaft import *
from feos.eos import *

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

fld = 'propane'
params = PcSaftParameters.from_json(['carbon dioxide', fld], 'gross2001.json')
kij = 0.116

params_kij = PcSaftParameters.new_binary(params.pure_records, kij)
pcsaft = EquationOfState.pcsaft(params_kij)
x1 = 0.5

phase_diagram = PhaseDiagram.bubble_point_line(pcsaft,np.array([x1, 1-x1])*MOL, min_temperature=200*KELVIN, npoints=501)
df = pd.DataFrame(phase_diagram.to_dict())
plt.plot(df['density liquid'],df.pressure/1e6,'b-')

phase_diagram = PhaseDiagram.dew_point_line(pcsaft,np.array([x1, 1-x1])*MOL, min_temperature=200*KELVIN, npoints=501)
df = pd.DataFrame(phase_diagram.to_dict())
plt.plot(df['density vapor'],df.pressure/1e6,'r-')

plt.ylabel(r'$p$ (MPa)')
plt.xlabel(r'$\rho$ ($\mathregular{mol/m^3}$)')

[g-bauer]

Yeah, that looks good. If you need specific pressures (and not a discretization in temperature), you can always generate the states by hand as I noted above.

Here is an example (using your code) that shows that both approaches lead to the same results:

from feos.si import *
from feos.pcsaft import *
from feos.eos import *

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

fld = 'propane'
params = PcSaftParameters.from_json(['carbon dioxide', fld], '../parameters/pcsaft/gross2001.json')
kij = 0.116

params_kij = PcSaftParameters.new_binary(params.pure_records, kij)
pcsaft = EquationOfState.pcsaft(params_kij)

x1 = 0.5
x = np.array([x1, 1-x1])

# Phase diagram for given pressures (between Pmin and 0.99*Pc)
pc = State.critical_point(pcsaft, x*MOL).pressure() * 0.99
pressures = SIArray1.linspace(0.1 * MEGA * PASCAL, pc, 20)
states = []
for pressure in pressures:
    if len(states) != 0:
        temperature = states[-1].liquid.temperature
        vapor_molefracs = states[-1].vapor.molefracs
    else:
        temperature = 200*KELVIN
        vapor_molefracs = None
    states.append(PhaseEquilibrium.bubble_point(
        pcsaft, 
        pressure, 
        liquid_molefracs=x, 
        vapor_molefracs=vapor_molefracs, 
        tp_init=temperature
    ))
densities = [state.liquid.density / (MOL / METER**3) for state in states]

# Phase diagram for given temperatures (between Tmin and Tc)
phase_diagram = PhaseDiagram.bubble_point_line(pcsaft,np.array([x1, 1-x1])*MOL, min_temperature=200*KELVIN, npoints=501)
df = pd.DataFrame(phase_diagram.to_dict())
plt.plot(df['density liquid'],df.pressure/1e6,'b-')

phase_diagram = PhaseDiagram.dew_point_line(pcsaft,np.array([x1, 1-x1])*MOL, min_temperature=200*KELVIN, npoints=501)
df = pd.DataFrame(phase_diagram.to_dict())
plt.plot(df['density vapor'],df.pressure/1e6,'r-')

plt.plot(densities, pressures / (MEGA * PASCAL), 'gx')

plt.ylabel(r'$p$ (MPa)')
plt.xlabel(r'$\rho$ ($\mathregular{mol/m^3}$)')

[g-bauer]

@RayleighSun is this issue resolved?

[RayleighSun]

@RayleighSun is this issue resolved?

yes, I will close it. Thanks!