Ctrl+K

Helmholtz EoS Parameter and Expression files

Helmholtz EoS Parameter and Expression files#

Helmholtz equations of state are difficult to implement in an equation oriented framework. Switching from the native temperature and density state variables to a more useful set and solving the phase equilibrium problem are difficult and may have multiple solutions, only one of which is physically meaningful. External functions allow for exact first and second derivatives to be calculated and used by the solver, while still allowing procedural code to be used for the solution. The external function also serve as decomposition making the models more robust by splitting off difficult-to-solve sub-problems to be solved separately.

To accommodate addition of substances using various forms of Helmholtz free energy models, Helmholtz free energy and transport properties are read into the external function library as AMPL NL-files. Basic parameters are read as a json file.

Generating the expression and parameter files needed is done via utilities provided by IDAES. The parameter and expression files originate from a Python script and an optional json master parameter file.

This section describes how to add or modify substances for the Helmholtz equation of state library.

Class#

The WriteParameters class creates the expression and parameter files needed by the external functions to define a substance. Parameters and expression are defined in a script and an optional main parameter file. Either custom or predefined expressions can be used.

class idaes.models.properties.general_helmholtz.helmholtz_parameters.WriteParameters(parameters)[source]#

This class allows users to set parameters and expressions to generate parameter and NL expression files that define a Helmholtz equation of state and optionally transport properties for a substance.

add(expressions)[source]#

This adds expressions to the to the object. These expressions are written to NL files to be used by external functions.

Parameters:

expressions (dict) – Dictionary where the key is an expression name and the value is a Pyomo expression.

Returns:

None

approx_sat_curves(trange)[source]#

_log.infos a table to verify that the approximate saturated density curves are correct. Since the approximate curves are used as an initial guess to the phase equilibrium problems, they are not directly testable.

Parameters:

trange (iterable) – temperature points in K

Returns:

list of liquid densities and list of vapor densities

Return type:

tuple

calculate_enthalpy(rho, T)[source]#

From the expressions provided, calculate enthalpy from density and temperature. This can be used for testing.

Parameters:

  • rho (float) – density in kg/m3

  • T (float) – temperature in K

Returns:

enthalpy in kJ/kg

Return type:

float

calculate_entropy(rho, T)[source]#

From the expressions provided, calculate entropy from density and temperature. This can be used for testing.

Parameters:

  • rho (float) – density in kg/m3

  • T (float) – temperature in K

Returns:

entropy in kJ/kg/K

Return type:

float

calculate_pressure(rho, T)[source]#

From the expressions provided, calculate pressure from density and temperature. This can be used for testing and calculating the critical pressure based on the critical temperature and density.

Parameters:

  • rho (float) – density in kg/m3

  • T (float) – temperature in K

Returns:

pressure in kPa

Return type:

float

calculate_reference_offset(delta, tau, s0, h0)[source]#

Given delta and tau for a reference state and the reference entropy and enthalpy, calculate the reference state offset parameters. Since temperature and density are independent of reference state, they can be calculated at a new reference state.

Parameters:

  • delta (float) – reduced density at new reference state

  • tau (float) – 1/reduced temperature at new reference state

  • s0 (float) – Entropy at new reference state [kJ/kg/K]

  • h0 (float) – Enthalpy at new reference state [kJ/kg]

Returns:

Offset for given reference state

Return type:

(tuple)

make_model(*args)[source]#

Make a Pyomo model used to define expression NL files. The arguments are strings for variables to create. The basic parameters are also added as parameters in the model.

Parameters:

(str) – Variables to add to the model

Returns:

Pyomo model with variables from args

Return type:

ConcreteModel

write(dry_run=False)[source]#

Write the parameter and expression files, and _log.info some diagnostics.

write_model(model, model_name, expressions=None)[source]#

Write an NL file and create an expression and variable map for the EoS model or just a variable map for transport property expressions where there is only one expression per file.

Parameters:

  • model (ConcreteModel) – a Pyomo model

  • model_name (str) – the model name used in the NL file

Returns:

NL file, expression map, variable map for EOS

Return type:

tuple

Master Parameter File#

While it is possible to the define equation of state and transport property expressions and parameters entirely in a Python script we will focus here on the use of a main parameter json file, which is read and converted into expressions and parameters that can be used by the external functions. The overall structure of the main json parameter file is given below. The details of each section will be filled as we progress through this documentation:

{
    "comp": "co2",
    "basic": {
        ...
    },
    "eos": {
        "reference": [
            "Span, R., and W. Wagner (1996). A New Equation of State for Carbon Dioxide",
            "    Covering the Fluid Region from the Triple-Point Temperature to 1100 K as",
            "    Pressures up to 800 MPa. Journal of Physical and Chemical Reference Data,",
            "    25, 1509."
        ],
        ...
    },
    "aux": {
        "reference": [
            "Span, R., and W. Wagner (1996). A New Equation of State for Carbon Dioxide",
            "    Covering the Fluid Region from the Triple-Point Temperature to 1100 K as",
            "    Pressures up to 800 MPa. Journal of Physical and Chemical Reference Data,",
            "    25, 1509."
        ],
        ...
    },
    "transport": {
        "thermal_conductivity": {
            "reference": [
                "Vesovic, V., W.A. Wakeham, G.A. Olchowy, J.V. Sengers, J.T.R. Watson, J.",
                "    Millat, (1990). The transport properties of carbon dioxide. J. Phys.",
                "    Chem. Ref. Data, 19, 763-808."
            ]
            ...
        },
        "viscosity": {
            "reference": [
                "Fenghour, A., W.A. Wakeham, V. Vesovic, (1998). The Viscosity of Carbon",
                "     Dioxide. J. Phys. Chem. Ref. Data, 27, 31-44."
            ]
            ...
        },
        "surface_tension": {
            "reference": [
                "Mulero, A., I. Cachadina, Parra, M., Recommended Correlations for the Surface ",
                "    Tension of Common Fluids, J. Phys. Chem. Ref. Data 41, 043105, (2012)."
            ],
            ...
        }
    }
}

By convention the main parameter file is named {comp}.json, where {comp} is the component name (e.g. co2.json); although, the main parameter file name can be anything.

Script#

An example of the most basic script for reading a main parameter file and generating expression and parameter files is given below. Predefined transport expressions are not yet available, so the script will also commonly define transport property expressions, and that will be covered later:

from idaes.models.properties.general_helmholtz.helmholtz_parameters import (
WriteParameters,
)


def main():
    """Generate parameter and expression files."""
    we = WriteParameters(parameters="r227ea.json")
    we.write()


if __name__ == "__main__":
    main()

Basic Parameters#

The basic parameters required are listed in the table below.

Parameter

Description

Units

R

specific gas constant (ideal gas constant/molecular weight)

kJ/kg/K

MW

molecular weight

g/mol

T_star

reducing temperature (usually critical temperature)

K

rho_star

reducing density (usually critical density)

kg/m3

Tc

critical temperature

K

rhoc

critical density

kg/m3

Pc

critical pressure (recalculated to match Tc and rhoc)

kPa

Tt

triple point temperature (can be approximate)

K

Pt

triple point pressure (can be approximate)

kPa

rhot_l

triple point liquid density (can be approximate)

kg/m3

rhot_v

triple point vapor density (can be approximate)

kg/m3

P_min

minimum pressure

kPa

P_max

maximum pressure

kPa

rho_max

maximum density

kg/m3

T_min

minimum temperature

K

T_max

maximum temperature

K

A truncated example of the main json parameter file is given below to show the form for basic parameters:

{
    "comp": "co2",
    "basic": {
        "R": 0.1889241,
        "MW": 44.0098,
        "T_star": 304.1282,
        "rho_star": 467.6,
        "Tc": 304.1282,
        "rhoc": 467.6,
        "Pc": 7377.3,
        "Tt": 216.592,
        "Pt": 517.95,
        "rhot_l": 1178.46,
        "rhot_v": 13.761,
        "P_min": 1e-9,
        "P_max": 5e5,
        "rho_max": 1600.0,
        "T_min": 216,
        "T_max": 1400
    },
    ...
}

Reference State#

Some properties such as enthalpy, entropy, and internal energy don’t have a well defined absolute value, and can only really be measured relative to a reference state. For the purpose of the parameter files, equations of state are initially implemented using whatever reference state is used in the original paper. An offset can be set in the parameter files to adjust the reference state if desired.

The calculate_reference_offset(delta, tau, s0, h0) method in the WriteParameters class can be used to calculate a new reference state. The temperature, density, reference enthalpy, and reference entropy are required. The temperature and density can be calculated using the IDAES property functions possibly using the expression writer class. The reference enthalpy and entropy are chosen arbitrarily. The reference state offset is dimensionless and is added to the \(n^\circ_1\) and \(n^\circ_2\) parameters in the ideal part for the dimensionless Helmholtz free energy corelation.

The following example shows how to set the reference state in the main parameter file:

{
    "comp": "co2",
    "basic": {
        ...
    },
    "eos": {
        ...
        "reference_state_offset": [
            -14.4979156224319,
            8.82013935801453
        ],
        ...
    },
    ...
}

Predefined Expressions#

Commonly used expressions are predefined and can be used by specifying expression types in the main parameter file.

Ideal Part of Dimensionless Helmholtz Free Energy#

Type 1

The predefined type 1 form of the ideal portion of the dimensionless Helmholtz free energy is shown below. Parameters are provided in the "eos" section of the main parameter file. This form of the expression requires a dictionary of the \(n^\circ_i$`\) parameters as n0 and the math gamma^circ$ parameters as g0. The last term index in the sum (\(h\)) is provided as last_term_ideal.

\[\log_e \delta + n^\circ_1 + n^\circ_2 \tau + n^\circ_3 \log_e \tau + \sum_{i = 4}^h n^\circ_i \log_e \left[ 1 - \exp(-\gamma^\circ_i \tau)\right]\]

A truncated example of a main parameter file is provided below, which uses this the type 1 form of the dimensionless Helmholtz free energy:

{
    "comp": "co2",
    "basic": {
        ...
    },
    "eos": {
        ...
        "n0": {
            "1": 8.37304456,
            "2": -3.70454304,
            "3": 2.5,
            "4": 1.99427042,
            "5": 0.62105248,
            "6": 0.41195293,
            "7": 1.04028922,
            "8": 0.08327678
        },
        "g0": {
            "4": 3.15163,
            "5": 6.11190,
            "6": 6.77708,
            "7": 11.32384,
            "8": 27.08792
        },
        "last_term_ideal": 8,
        "phi_ideal_type": 1,
        ...
    },
    ...
}

Type 2

Parameters for the type 2 version are similar to type 1, but the last_term_ideal parameter is a list [h1, h2].

\[\log_e \delta + n^\circ_1 + n^\circ_2 \tau + n^\circ_3 \log_e \tau + \sum_{i = 4}^{h_1} n^\circ_i \tau^{\gamma^\circ_i} + \sum_{i = h_1 + 1}^{h_2} n^\circ_i \log_e \left[ 1 - \exp(-\gamma^\circ_i \tau)\right]\]

Type 3

Parameters for the type 3 version are similar to type 1.

\[\log_e \delta + n^\circ_1 + n^\circ_2 \tau + n^\circ_3 \log_e \tau + \sum_{i = 4}^{h} n^\circ_i \tau^{\gamma^\circ_i}\]

Residual Part of Dimensionless Helmholtz Free Energy#

Predefined residual term expressions are discussed here. The last_term_residual parameter is a list [h1, h2, …]. Equation parameters are given by dictionaries in the "eos" section. Greek letter parameters are given by the parameters \(\alpha\) = a, \(\beta\) = b, \(\varepsilon\) = e, and \(\gamma\) = g.

The example below shows how to use a predefined residual dimensionless Helmholtz free energy expression:

{
    "comp": "co2",
    ...
    "eos": {
        ...
        "c": {
            "8": 1,
            ...
            "34": 6
        },
        "d": {
            "1": 1,
            ...
            "39": 3
        },
        "t": {
            "1": 0.00,
            ...
            "39": 3.00
        },
        "n": {
            "1": 3.88568232031610e-01,
            ...
            "42": 5.50686686128420e-02
        },
        "a": {
            "35": 25,
            ...
            "39": 20
        },
        "b": {
            "35": 325,
            ...
            "39": 275
        },
        "g": {
            "35": 1.16,
            ...
            "39": 1.22
        },
        "e": {
            "35": 1,
            "36": 1,
            "37": 1,
            "38": 1,
            "39": 1
        },
        "last_term_residual": [7, 34, 39],
        "phi_residual_type": 2
        ...
    }
    ...
}

Type 1

\[\phi^r(\delta, \tau) = \sum^{h_1}_1 n_i \delta^{d_i} \tau^{t_i} + \sum_{h_1 + 1}^{h_2} n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{c_i})\]

Type 2

\[\phi^r(\delta, \tau) = \sum^{h_1}_{i=1} n_i \delta^{d_i} \tau^{t_1} + \sum_{i =h_1 + 1}^{h_2} n_i \delta^{d_i} \tau^{t_1} \exp(-\delta^{c_i}) + \sum_{i = h_2 + 1}^{h_3} n_i \delta^{d_i} \tau^{t_1} \exp\left[-\alpha_i(\delta - \varepsilon_i)^2 - \beta_i(\tau - \gamma_i)^2 \right]\]

Type 3

\[\phi^r(\delta, \tau) = \sum^{h_1}_{i=1} n_i \delta^{d_i} \tau^{t_1} + \sum_{i=h_1 + 1}^{h_2} n_i \delta^{d_i} \tau^{t_1} \exp(-\delta^{c_i}) \sum_{i=h_2 + 1}^{h_3} n_i \delta^{d_i} \tau^{t_1} \exp(-\delta^{c_i}) \exp(-\tau^{b_i})\]

Type 4

\[\phi^r(\delta, \tau) = \sum^{h_0}_{i=1} n_i \delta^{d_i} \tau^{t_1} + \sum^m_{j=1} \left[ \exp(-\delta^j) \sum^{h_j}_{i = h_{j-1} + 1} n_i \delta^{d_i} \tau^{t_1} \right]\]

Approximate Saturated Reduced Density#

To assist in the equilibrium calculation, curves for approximate saturated liquid and vapor density are required. These curves are typically provided in the papers where the equations of state are defined. There are two common forms for these curves, which we call type 1 and type 2 for convenience.

An example main parameter file entry for the approximate density curves is given below:

{
    "comp": "co2",
    ...
    "aux": {
        "reference": [
            "Span, R., and W. Wagner (1996). A New Equation of State for Carbon Dioxide",
            "    Covering the Fluid Region from the Triple-Point Temperature to 1100 K as",
            "    Pressures up to 800 MPa. Journal of Physical and Chemical Reference Data,",
            "    25, 1509."
        ],
        "delta_l_sat_approx": {
            "c": 1,
            "n": {
                "1": 1.9245108,
                "2": -0.62385555,
                "3": -0.32731127,
                "4": 0.39245142
            },
            "t": {
                "1": 0.34,
                "2": 0.5,
                "3": 1.6666666666666667,
                "4": 1.8333333333333333
            },
            "type": 2
        },
        "delta_v_sat_approx": {
            "c": 1,
            "n": {
                "1": -1.7074879,
                "2": -0.82274670,
                "3": -4.6008549,
                "4": -10.111178,
                "5": -29.742252
            },
            "t": {
                "1": 0.340,
                "2": 0.5,
                "3": 1.0,
                "4": 2.3333333333333333,
                "5": 4.666666666666667
            },
            "type": 2
        }
    },
    ...
}

Type 1

\[\delta = c + \sum_{i=1}^h n_i \left(1 - \frac{T}{T_c} \right)^{t_i}\]

Type 2

\[\delta = c \exp \left[ \sum_{i=0}^h \left(1 - \frac{T}{T_c} \right)^{t_i} \right]\]

Surface Tension#

The surface tension external functions provide surface tension in units of mN/m. There is only one form for surface tension and it is only a function of temperature. The parameters required are Tc [K], s [mN/m] and n [dimensionless]. Tc is provided again in case the critical temperature differs slightly from the equation of state parameter. The number of terms is inferred from the s dictionary. This version of the surface tension expression maxes out at the critical surface tension, so as not to cause numerical issues at temperatures above critical.

\[\sigma = \sum_{i = 0}^h s_i \max\left(1 - \frac{T}{T_c}, 0 \right)^{n_i}\]

An example surface tension section is given below:

{
    "comp": "co2",
    "transport": {
        ...
        "surface_tension": {
            "reference": [
                "Mulero, A., I. Cachadina, Parra, M., Recommended Correlations for the Surface ",
                "    Tension of Common Fluids, J. Phys. Chem. Ref. Data 41, 043105, (2012)."
            ],
            "Tc": 304.1282,
            "s": {
                "0": 78.63
            },
            "n": {
                "0": 2.471
            },
            "type": 1
        }
    }
}

Viscosity#

Currently, there are no predefined viscosity expressions.

Thermal Conductivity#

Currently, there are no predefined viscosity expressions.

Custom Expressions#

Although any custom expressions can be defined rather than using the predefined ones, we will focus here on viscosity and thermal conductivity since predefined expressions are not available from them yet. The water example demonstrates all the concepts.

Variables#

Custom expressions are written on Pyomo models. With the delta (reduced density) and tau (1/reduced temperature). The easiest way to add an expression is to write a rule for it then use the add() method to add the expression to the WriteParameters class.

Units#

Most expression used are dimensionless. Only viscosity (microPa*s), thermal conductivity (mW/m/K) and surface tension (mN/m) have units.

Property Functions#

If expressions require calculation of other properties (e.g. pressure or heat capacity) they can call other property functions, as demonstrated in the example up next.

Example#

The example for water below shows how custom expressions can be defined. Note that heat capacity, viscosity, and isothermal compressibility are used in the thermal conductivity expression by calling property functions:

import math
import pyomo.environ as pyo
from pyomo.common.fileutils import find_library
from idaes.core.util.math import smooth_max
from idaes.models.properties.general_helmholtz.helmholtz_parameters import (
    WriteParameters,
)

def thermal_conductivity_rule(m):
    """Thermal conductivity rule

    International Association for the Properties of Water and Steam (2011).
        IAPWS R15-11, "Release on the IAPWS Formulation 2011 for the
        Thermal Conductivity of Ordinary Water Substance,"
        URL: http://iapws.org/relguide/ThCond.pdf
    """

    L0 = {
        0: 2.443221e-3,
        1: 1.323095e-2,
        2: 6.770357e-3,
        3: -3.454586e-3,
        4: 4.096266e-4,
    }

    L1 = {
        (0, 0): 1.60397357,
        (1, 0): 2.33771842,
        (2, 0): 2.19650529,
        (3, 0): -1.21051378,
        (4, 0): -2.7203370,
        (0, 1): -0.646013523,
        (1, 1): -2.78843778,
        (2, 1): -4.54580785,
        (3, 1): 1.60812989,
        (4, 1): 4.57586331,
        (0, 2): 0.111443906,
        (1, 2): 1.53616167,
        (2, 2): 3.55777244,
        (3, 2): -0.621178141,
        (4, 2): -3.18369245,
        (0, 3): 0.102997357,
        (1, 3): -0.463045512,
        (2, 3): -1.40944978,
        (3, 3): 0.0716373224,
        (4, 3): 1.1168348,
        (0, 4): -0.0504123634,
        (1, 4): 0.0832827019,
        (2, 4): 0.275418278,
        (3, 4): 0.0,
        (4, 4): -0.19268305,
        (0, 5): 0.00609859258,
        (1, 5): -0.00719201245,
        (2, 5): -0.0205938816,
        (3, 5): 0.0,
        (4, 5): 0.012913842,
    }
    delta = m.delta
    tau = m.tau
    big_lam = 177.8514
    qdb = 1.0 / 0.40
    nu = 0.630
    gamma = 1.239
    xi0 = 0.13
    big_gam0 = 0.06
    Tbr = 1.5
    Tb = 1 / tau
    lib_path = find_library("general_helmholtz_external")
    m.cp = pyo.ExternalFunction(library=lib_path, function="cp")
    m.cv = pyo.ExternalFunction(library=lib_path, function="cv")
    m.mu = pyo.ExternalFunction(library=lib_path, function="mu")
    m.itc = pyo.ExternalFunction(
        library=lib_path, function="itc"
    )  # isothermal compressibility [1/MPa]
    deltchi = smooth_max(
        delta**2
        * (m.itc("h2o", delta, tau) - m.itc("h2o", delta, 1 / Tbr) * Tbr / Tb)
        * m.Pc
        / 1000,
        0,
        eps=1e-8,
    )
    xi = xi0 * (deltchi / big_gam0) ** (nu / gamma)
    y = qdb * xi
    kappa = m.cp("h2o", delta, tau) / m.cv("h2o", delta, tau)
    Z = (
        2.0
        / math.pi
        / y
        * (
            ((1 - 1.0 / kappa) * pyo.atan(y) + 1.0 / kappa * y)
            - (1 - pyo.exp(-1 / (1.0 / y + y**2 / 3 / delta**2)))
        )
    )
    cpb = m.cp("h2o", delta, tau) / m.R
    mub = m.mu("h2o", delta, tau)

    lam2 = big_lam * delta * Tb * cpb / mub * Z
    return lam2 + pyo.sqrt(1.0 / m.tau) / sum(
        L0val * m.tau**i for i, L0val in L0.items()
    ) * pyo.exp(
        m.delta
        * sum(
            (m.tau - 1) ** i * sum(L1[i, j] * (m.delta - 1) ** j for j in range(0, 6))
            for i in range(0, 5)
        )
    )


def viscosity_rule(mvisc):
    """Viscosity calculation for water

    International Association for the Properties of Water and Steam (2008).
        IAPWS R12-08, "Release on the IAPWS Formulation 2008 for the Viscosity
        of Ordinary Water Substance, URL: http://iapws.org/relguide/visc.pdf
    """

    H0 = {0: 1.67752, 1: 2.20462, 2: 0.6366564, 3: -0.241605}
    H1 = {
        (0, 0): 5.20094e-1,
        (1, 0): 8.50895e-2,
        (2, 0): -1.08374,
        (3, 0): -2.89555e-1,
        (4, 0): 0.0,
        (5, 0): 0.0,
        (0, 1): 2.22531e-1,
        (1, 1): 9.99115e-1,
        (2, 1): 1.88797,
        (3, 1): 1.26613,
        (4, 1): 0.0,
        (5, 1): 1.20573e-1,
        (0, 2): -2.81378e-1,
        (1, 2): -9.06851e-1,
        (2, 2): -7.72479e-1,
        (3, 2): -4.89837e-1,
        (4, 2): -2.57040e-1,
        (5, 2): 0.0,
        (0, 3): 1.61913e-1,
        (1, 3): 2.57399e-1,
        (2, 3): 0.0,
        (3, 3): 0.0,
        (4, 3): 0.0,
        (5, 3): 0.0,
        (0, 4): -3.25372e-2,
        (1, 4): 0.0,
        (2, 4): 0.0,
        (3, 4): 6.98452e-2,
        (4, 4): 0.0,
        (5, 4): 0.0,
        (0, 5): 0.0,
        (1, 5): 0.0,
        (2, 5): 0.0,
        (3, 5): 0.0,
        (4, 5): 8.72102e-3,
        (5, 5): 0.0,
        (0, 6): 0.0,
        (1, 6): 0.0,
        (2, 6): 0.0,
        (3, 6): -4.35673e-3,
        (4, 6): 0.0,
        (5, 6): -5.93264e-4,
    }

    return (
        1e2
        * pyo.sqrt(1.0 / mvisc.tau)
        / sum(H0val * mvisc.tau**i for i, H0val in H0.items())
        * pyo.exp(
            mvisc.delta
            * sum(
                (mvisc.tau - 1) ** i
                * sum(H1[i, j] * (mvisc.delta - 1) ** j for j in range(0, 7))
                for i in range(0, 6)
            )
        )
    )


def main():
    """Generate property and expression files."""

    we = WriteParameters(parameters="h2o.json")
    we.add(
        {
            "viscosity": viscosity_rule,
            "thermal_conductivity": thermal_conductivity_rule,
        }
    )
    we.write()


if __name__ == "__main__":
    main()
Contents