Error when computing vle of polyethylene + supercritical ethylene

[Paiding]

I tried to compute the phase equilibrium of PE and ethylene with PC-SAFT EOS to reproduce the results of [T. Spyriouni, I.G. Economou / Polymer 46 (2005) 10772–10781]

I modified the code in examples/pcsaft/PhaseDiagrmBinary.ipynb for supercritical mixture and is like this:

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

import matplotlib.pyplot as plt
import numpy as np

ethylene = PureRecord(
    identifier=Identifier(
        cas='74-85-1',
        name='ethylene',
        iupac_name='ethene',
        smiles='C=C',
        inchi='InChI=1/C2H4/c1-2/h1-2H2',
        formula='C2H4'
    ),
    molarweight=28.05,
    model_record=PcSaftRecord(
        m=1.593,
        sigma=3.445,
        epsilon_k=176.47,
    )
)

polyethylene = PureRecord(
    identifier=Identifier(
        cas='9002-88-4"',
        name='polyethylene',
        iupac_name='whatever',
        smiles='whatever',
        inchi='Whatever',
        formula='[C2H4]n'
    ),
    molarweight=99000,
    model_record=PcSaftRecord(
        m=2603.7,
        sigma=3.75,
        epsilon_k=252.00,
    )
)

params = PcSaftParameters.new_binary([ethylene, polyethylene], binary_record=0.0427)

saft = EquationOfState.pcsaft(params)
dia_p = PhaseDiagram.binary_vle(saft, 423.15*KELVIN)

But I got a error that

dia_p = PhaseDiagram.binary_vle(saft, 423.15*KELVIN)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: System is supercritical.

I also tried the mixture of supercritical ethylene with shorter alkanes, FeOs pcsaft works well. Is it a error solving equations for polymer’s properties?

[prehner]

Thank you for bringing this to our attention and the comprehensive example! While PC-SAFT can and has been used to describe properties of polymers, it was not a high priority of ours during the development of FeOs. The special circumstances like the high pressures and the extremely low (molar) concentrations of the polymer pose challenges on standard phase equilibrium calculations. The calculation of properties (e.g. activity or fugacity coefficients) can be expected to work perfectly fine for polymers. However, you would need a dedicated algorithm to determine the cloud or shadow points of polymer/monomer mixtures.

[Paiding]

I'm very grateful for your helpful reply. My tutor suggested me to draw a μ-x curve for each NTP states to find the phase separation point. However, when I tried to get the chemical_potential value from one state:

state_npt = State(
    saft,
    temperature=temp*KELVIN,
    pressure=press*BAR,
    molefracs=np.array([0.5,0.5]),
    total_moles=MOL,
    density_initialization='liquid'
)
state_npt.chemical_potential()

It retruned pyo3_runtime.PanicException: No ideal gas model initialized! in my system. (for simpler systems, too)
It seems that the expressions about chemical potential in tutorials (pcsaft_state.ipynb/core_user_defined_eos.ipynb) are not suitable for current 0.6.1 FeOs. Would you please show us where to set the 'ideal gas model' or how to get chemical potential for a state?

[g-bauer]

Hello Paiding,

sorry for the poor documentation, it's something we are working on.

First of all, you don't need an ideal gas model if you are only interested in the residual chemical potential. You can use state_npt.chemical_potential(Contributions.Residual) to achieve that.

If you want to calculate the total chemical potential (the default when no Contribution object is used as argument), you have three options at the moment.

1. define the model for the de Broglie wavelength yourself in Python
2. use a DIPPR correlation (we provide a Dippr object, but we don't have parameters for the substances you want to use)
3. use a group-contribution method (we provide Joback for that)

Here is an example of manually building the group contributions. Please check if those values make sense, it's merely for demonstration purposes (i.e. I just use CH2 records for polyethylene).

from feos.si import *
from feos.pcsaft import *
from feos.eos import *
from feos.joback import Joback, SegmentRecord, JobackRecord

import matplotlib.pyplot as plt
import numpy as np

ethylene = PureRecord(
    identifier=Identifier(name='ethylene'),
    molarweight=28.05,
    model_record=PcSaftRecord(
        m=1.593,
        sigma=3.445,
        epsilon_k=176.47,
    )
)

polyethylene = PureRecord(
    identifier=Identifier(name='polyethylene'),
    molarweight=99000,
    model_record=PcSaftRecord(
        m=2603.7,
        sigma=3.75,
        epsilon_k=252.00,
    )
)

# Build Joback segments + parameters
joback_ethylene = ChemicalRecord(
    identifier=Identifier(name='ehtylene'),
    segments=['=CH2']*2 # list of segments that make up the molecule
)

joback_polyethylene = ChemicalRecord(
    identifier=Identifier(name='polyehtylene'),
    segments=['CH2']*int(99000 // 14.02658)
)

# Taken from parameters/ideal_gas/joback1987.json
dch2r = SegmentRecord(
    '=CH2', 
    molarweight=14.02658, 
    model_record=JobackRecord(23.6, -0.0381, 0.000172, -1.03e-07, 0.0)
)

ch2r = SegmentRecord(
    'CH2', 
    molarweight=14.02658, 
    model_record=JobackRecord(-0.909, 0.095, -5.44e-05, 1.19e-08, 0.0)
)

joback_parameters = Joback.from_segments(
    chemical_records=[joback_ethylene, joback_polyethylene], 
    segment_records=[dch2r, ch2r]
)

saft_parameters = PcSaftParameters.new_binary([ethylene, polyethylene], binary_record=0.0427)
# build eos including ideal gas model
saft = EquationOfState.pcsaft(saft_parameters).joback(joback_parameters)

state_npt = State(
    saft,
    temperature=300*KELVIN,
    pressure=BAR,
    molefracs=np.array([0.5, 0.5]),
    density_initialization='liquid'
)

print("mu_res", state_npt.chemical_potential(Contributions.Residual))
print("mu_ig ", state_npt.chemical_potential(Contributions.IdealGas))
print("mu_tot", state_npt.chemical_potential()) # equivalent to state_npt.chemical_potential(Contributions.Total)

Building parameters for GC methods is discussed in detail here. Let me know if you have any questions.