Incorrect densities for a Peng-Robinson fluid when it should be in the liquid phase

[miloknowles]

Problem description

When creating a Solution for methanol using a Peng-Robinson equation of state, the incorrect density (off by orders of magnitude) is returned at temperatures and pressure where we'd expect methanol to be in a liquid phase. For example, at a TP of (300, 101325), I see a density of ~1.3 kg/m3, whereas liquid methanol would be around ~780 kg/m3.

Digging into the C++ code, it appears that there is some logic to "force" a certain solution to the cubic equation of state when multiple are found. However, I don't see any Python interface for doing that. In general, it's unclear whether a solution using the Peng-Robinson equation of state will reliable detect whether a fluid is in the liquid or vapor phase, and there doesn't seem to be a way to "help" the solver find the right density. If you can point me in the right direction, I'd be happy to work on a PR that addresses this issue.

Steps to reproduce

Using the following YAML file:

# methanol.yaml
generator: cti2yaml
cantera-version: 2.5.0

phases:
- name: methanol
  thermo: Peng-Robinson
  elements: [C, H, O]
  species: [CH3OH]
  state:
    T: 200
    P: 1.01325e+05

species:
- name: CH3OH
  composition: {C: 1, H: 4, O: 1}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 1000.0, 3500.0]
    data:
    - [5.71539582, -0.0152309129, 6.52441155e-05, -7.10806889e-08, 2.61352698e-11,
      -2.56427656e+04, -1.50409823]
    - [1.78970791, 0.0140938292, -6.36500835e-06, 1.38171085e-09, -1.1706022e-13,
      -2.53748747e+04, 14.5023623]
    note: L8/88
  equation-of-state:
    model: Peng-Robinson
    a: 1.0138872136786787
    b: 4.041811555663748e-05
    acentric-factor: 0.5625

I load a Solution using Python:

spec = ct.Solution('methanol.yaml')

# Boiling point, used below
Tb = 64.96 + 273.15

# Source of comparison:
# https://www.engineeringtoolbox.com/methanol-density-specific-weight-temperature-pressure-d_2091.html

# Expect to be a liquid:
spec.TP = 300, 101325
print("Density @ 300K:", spec.density_mass) # Expected: 784.3 kg/m3

# Expect to be a liquid:
spec.TP = Tb - 10, 101325
print("Density @ Tb-10:", spec.density_mass)

# Expect to be a gas:
spec.TP = 400, 101325
print("Density @ 400K:", spec.density_mass) # Expected: 0.976 kg/m3

spec.TP = 500, 101325
print("Density @ 500K:", spec.density_mass) # Expected: 0.774 kg/m3

Printout:

Density @ 300K: 1.3016220571692458
Density @ Tb-10: 1.190107980566739
Density @ 400K: 0.9762144154660121
Density @ 500K: 0.7809704748372757

Behavior

As you can see, the density is very low (indicating a gas) even at temperatures below the boiling point of methanol. Expected values are shown as comments to the right of each line. I haven't debugged the C++ code to confirm this, but I suspect that there are multiple roots of the cubic polynomial, and the one corresponding to a gas phase is being returned.

System information

• Cantera 3.0.0, installed with poetry
• Mac M2 OSX
• Python 3.12.3

Attachments

Additional context

[wandadars]

In my experience, the Peng Robinson equation of state is mostly meant for gaseous solutions. The MixtureFugacityTP class is the one responsible for deciding phase information. Unfortunately, I believe that class is incomplete and mostly the gas information is what it uses.

[wandadars]

I think a better handling of the detection of the saturation temperature is what would be needed. The Gibbs balance between the liquid and vapor phases would allow for a better selection of the vapor or liquid solution. One thing is that I don't think any two-phase solutions are permitted currently.

[ischoegl]

@wandadars is correct - in its current form, Cantera is not set up for two-phase Solutions (outside of pure substances). Peng-Robinson et al are currently implemented to capture real gas behavior only. There are, however, some known glitches with the cubic solver, see #1699 or #1157, so logic improvements would be welcome. Beyond, all phases are implemented in C++ only; Python only serves as a more user-friendly interface, and does not provide access to all internal functions.

PS: apart from that, it does appear that some exception handling would be beneficial to indicate when the state can no longer be represented correctly.