Cubic Equations of State

Introduction

CoolProp (as of version 6) comes with two standard cubic equations of state: Soave-Redlich-Kwong (SRK) and Peng-Robinson (PR). These two equations of state can be expressed in a common, generalized form:

\[p = \frac{RT}{v-b} + \frac{a}{(v+\Delta_1b)(v+\Delta_2b)}\]

where for pure fluids, a and b are not composition dependent, whereas for mixtures, they have composition dependence.

The motivations for the use of cubic EOS are twofold:

  • The only required information for the EOS are \(T_c\), \(p_c\), and the acentric factor of the pure fluids
  • They are much more computationally efficient (see below)

Caveats

Warning

NOT ALL PROPERTIES ARE AVAILABLE AS INPUTS/OUTPUTS

Only a limited subset of properties are available currently. You can do:

  • Flash calculations with TP, PQ, DT, QT inputs
  • Calculation of mixture critical point(s)
  • Calculation of some mixture flashes

Pure Fluids

Usage

As a user, in the high-level interface, all that you have to do to evaluate properties using the SRK or PR backends is to change the backend name.

In [1]: import CoolProp.CoolProp as CP

# The multi-parameter Helmholtz backend
In [2]: CP.PropsSI("T","P",101325,"Q",0,"HEOS::Propane")
Out[2]: 231.03621464431768

# SRK
In [3]: CP.PropsSI("T","P",101325,"Q",0,"SRK::Propane")
Out[3]: 231.27041592485747

# Peng-Robinson
In [4]: CP.PropsSI("T","P",101325,"Q",0,"PR::Propane")
Out[4]: 230.95724200364225

The same holds for the low-level interface:

In [5]: import CoolProp.CoolProp as CP

In [6]: AS = CP.AbstractState("SRK", "Propane"); AS.update(CP.QT_INPUTS, 0, 300); print(AS.p())
1008752.02868

In [7]: AS = CP.AbstractState("PR", "Propane"); AS.update(CP.QT_INPUTS, 0, 300); print(AS.p())
997556.081642

All the fluids available in CoolProp are also available through the cubic equations of state. The fluids can be extended according to the analysis shown below

Speed

The increase in speed for evaluating properties from a cubic EOS is one of the primary motivations. Here we show an example of the speedup to VLE calculations:

In [8]: import CoolProp.CoolProp as CP

# The multi-parameter Helmholtz backend
In [9]: AS = CP.AbstractState("HEOS", "Propane")

In [10]: %timeit AS.update(CP.QT_INPUTS, 0, 300)
10000 loops, best of 3: 30.6 us per loop

# The cubic SRK backend
In [11]: AS = CP.AbstractState("SRK", "Propane")

In [12]: %timeit AS.update(CP.QT_INPUTS, 0, 300)
The slowest run took 5.29 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 4.91 us per loop

And here, we run the PT flash

In [13]: import CoolProp.CoolProp as CP

# The multi-parameter Helmholtz backend
In [14]: AS = CP.AbstractState("HEOS", "Propane")

In [15]: AS.specify_phase(CP.iphase_gas)

In [16]: %timeit AS.update(CP.PT_INPUTS, 101325, 300)
100000 loops, best of 3: 12.4 us per loop

# The cubic SRK backend
In [17]: AS = CP.AbstractState("SRK", "Propane")

In [18]: AS.specify_phase(CP.iphase_gas)

In [19]: %timeit AS.update(CP.PT_INPUTS, 101325, 300)
The slowest run took 13.82 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 1 us per loop

As you can see, the speed difference is quite significant

Mixtures

Interaction Parameters

For mixtures, cubic EOS (and their modifications) are heavily used. Cubic EOS allow for reasonable prediction of VLE with only one adjustable parameter, at least for reasonable mixtures. CoolProp does not come with any values for the \(k_{ij}\) parameter, but it is straightforward to add it yourself.

Warning

The ability to adjust interaction parameters is ONLY available in the low-level interface

Warning

When you call set_binary_interaction_double, it only applies to the given instance of the AbstractState

In [20]: import CoolProp.CoolProp as CP

In [21]: AS = CP.AbstractState("SRK", "Methane&Ethane")

In [22]: AS.set_mole_fractions([0.5, 0.5])

In [23]: AS.update(CP.QT_INPUTS, 0, 120); print(AS.p())
104467.462374

In [24]: AS.set_binary_interaction_double(0,1,"kij",-0.05)

In [25]: AS.update(CP.QT_INPUTS, 0, 120); print(AS.p())
82034.1034555

Critical Points

According to a forthcoming paper from Bell et al., it is possible to calculate critical points of mixtures from cubic EOS. For instance, here is how to calculate all the critical points (there can be more than one) that are found for an equimolar mixture of methane and ethane:

In [26]: import CoolProp.CoolProp as CP

In [27]: AS = CP.AbstractState("SRK", "Methane&Ethane")

In [28]: AS.set_mole_fractions([0.5, 0.5])

In [29]: pts = AS.all_critical_points()

In [30]: [(pt.T, pt.p, pt.rhomolar, pt.stable) for pt in pts]
Out[30]: [(266.3194710014874, 6859304.778393156, 7895.129640980038, True)]

Detailed Example

Here we plot phase envelopes and critical points for an equimolar methane/ethane mixture

(Source code, png, .pdf)

../_images/Cubics-1.png

Adding your own fluids

The cubic fluids in CoolProp are defined based on a JSON format, which could yield something like this for a FAKE (illustrative) fluid. For instance if we had the file fake_fluid.json (download it here: fake_fluid.json):

[
  {
    "CAS": "000-0-00", 
    "Tc": 400.0, 
    "Tc_units": "K", 
    "acentric": 0.1, 
    "aliases": [
    ], 
    "molemass": 0.04, 
    "molemass_units": "kg/mol", 
    "name": "FAKEFLUID", 
    "pc": 4000000.0, 
    "pc_units": "Pa"
  }
]

The JSON-formatted fluid information is validated against the JSON schema, which can be obtained by calling the get_global_param_string function:

In [31]: import CoolProp.CoolProp as CP

# Just the first bit of the schema (it's a bit large; see below for the complete schema)
In [32]: CP.get_global_param_string("cubic_fluids_schema")[0:60]
Out[32]: u'{\n  "$schema": "http://json-schema.org/draft-04/schema#",\n  '

A JSON schema defines the structure of the data, and places limits on the values, defines what values must be provided, etc.. There are a number of online tools that can be used to experiment with and validate JSON schemas: http://www.jsonschemavalidator.net , http://jsonschemalint.com/draft4/.

In [33]: import CoolProp.CoolProp as CP, json

# Normally you would read it from a file!
In [34]: fake_fluids = [
   ....:                   {
   ....:                     "CAS": "000-0-00",
   ....:                     "Tc": 400.0,
   ....:                     "Tc_units": "K",
   ....:                     "acentric": 0.1,
   ....:                     "aliases": [
   ....:                     ],
   ....:                     "molemass": 0.04,
   ....:                     "molemass_units": "kg/mol",
   ....:                     "name": "FAKEFLUID",
   ....:                     "pc": 4000000.0,
   ....:                     "pc_units": "Pa"
   ....:                   }
   ....:               ]
   ....: 

# Adds fake fluid for both SRK and PR backends since they need the same parameters from the same internal library
In [35]: CP.add_fluids_as_JSON("PR", json.dumps(fake_fluids))

In [36]: CP.PropsSI("T","P",101325,"Q",0,"SRK::FAKEFLUID")
Out[36]: 247.0950826605273

Complete fluid schema

For completeness, here is the whole JSON schema used for the cubic backends:

In [37]: import CoolProp.CoolProp as CP

In [38]: print(CP.get_global_param_string("cubic_fluids_schema"))
{
  "$schema": "http://json-schema.org/draft-04/schema#",
  "title": "Fluid for cubic EOS in CoolProp",
  "items": {
    "properties": {
      "name": {
        "description": "Name of the fluid",
        "type": "string"
      },
      "CAS": {
        "description": "CAS registry number of the fluid",
        "type": "string"
      },
      "Tc": {
        "description": "Critical temperature",
        "type": "number",
        "minimum": 0.1,
        "maximum": 20000
      },
      "pc": {
        "description": "Critical pressure",
        "type": "number",
        "minimum": 0,
        "maximum": 500000000
      },
      "acentric": {
        "description": "Acentric factor",
        "type": "number",
        "minimum": -10,
        "maximum": 10
      },
      "molemass": {
        "description": "Molar mass",
        "type": "number",
        "minimum": 0,
        "maximum": 1
      },
      "molemass_units": {
        "description": "Units of the molar mass provided",
        "enum": [
          "kg/mol"
        ]
      },
      "pc_units": {
        "description": "Units of the critical pressure",
        "enum": [
          "Pa"
        ]
      },
      "Tc_units": {
        "description": "Units of the critical temperature",
        "enum": [
          "K"
        ]
      },
      "aliases": {
        "type": "array",
        "items": {
          "type": "string"
        },
        "minItems": 0
      },
      "alpha": {
        "description": "The alpha function being used",
        "properties": {
          "type":{
            "enum": [
              "Twu",
              "Mathias-Copeman",
              "default"
            ]
          },
          "c":{
              "type": "array",
              "items": {
                "type": "number"
              },
              "minItems": 3,
              "maxItems": 3
          }
        }
      }
    },
    "required": [
      "name",
      "CAS",
      "Tc",
      "Tc_units",
      "pc",
      "pc_units",
      "acentric",
      "molemass",
      "molemass_units"
    ]
  },
  "type": "array"
}