#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2024 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
#################################################################################
"""This provides functions to create and check Helmholtz equation of state
parameter and expression files.
"""
import logging
import json
import numpy
import pyomo.environ as pyo
from idaes.models.properties.general_helmholtz.expressions import (
phi_residual_types,
phi_ideal_types,
delta_sat_types,
surface_tension_types,
)
_log = logging.getLogger("idaes.helmholtz_parameters")
def _parse_int_key(pairs):
"""Since json can only store dictionary keys as strings this converts the string
keys to integers or integer tuples if possible
"""
d = {}
for x in pairs:
try:
if x[0][0] == "(" and x[0][-1] == ")":
d[tuple(map(int, map(str.strip, x[0][1:-1].split(","))))] = x[1]
else:
d[int(x[0])] = x[1]
except ValueError:
d[x[0]] = x[1]
return d
[docs]
class WriteParameters(object):
"""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.
"""
# Dictionary used to track the location of variables in NL files.
variables = {
"delta": 0,
"tau": 1,
"p": 2,
"T": 3,
}
# Dictionary used to track the location of expressions in NL files.
expressions = {
# phi is dimensionless Helmholtz free energy
"phii": 0, # ideal part of phi(delta, tau)
"phii_d": 1, # ideal part of phi partial wrt delta
"phii_dd": 2, # ideal part of phi partial wrt delta and delta
"phii_t": 3, # ideal part of phi partial wrt tau
"phii_tt": 4, # ideal part of phi partial tau delta and tau
"phii_dt": 5, # ideal part of phi partial wrt delta and tau
"phir": 6, # residual part of phi(delta, tau)
"phir_d": 7, # residual part of phi partial wrt delta
"phir_dd": 8, # residual part of phi partial wrt delta and delta
"phir_t": 9, # residual part of phi partial wrt tau
"phir_tt": 10, # residual part of phi partial wrt tau and tau
"phir_dt": 11, # residual part of phi partial wrt delta and tau
# Initial guess on saturated liquid and vapor reduced density
"delta_v_sat_approx": 12, # approximate delta_sat_vapor(tau)
"delta_l_sat_approx": 13, # approximate delta_sat_liquid(tau)
}
# A dictionary of optional expressions for transport properties, the value
# in the dictionary is a short name used in the NL file names.
optional_expressions = {
"thermal_conductivity": "tcx",
"surface_tension": "st",
"viscosity": "visc",
}
def __init__(
self,
parameters,
):
"""Create a parameter writer object. Parameters can either be specified in
dictionary or read from a json file.
Args:
parameters (str|dict): Either a dictionary of parameters or a path to
a json file.
Returns:
WriteParameters: Object used to write Helmholtz EoS parameter and expression files
"""
# A list that will store the names of defined expressions
self.has_expression = []
# if parameters is a string read from json
if isinstance(parameters, str):
with open(parameters, "r") as fp:
parameters = json.load(fp, object_pairs_hook=_parse_int_key)
# we may need to access parameters that are provided but not stored
# directly as attributes, so store a link to the parameter dictionary
self.parameters = parameters
# Store the basic parameters as attributes.
self.comp = parameters["comp"]
self.R = parameters["basic"]["R"]
self.MW = parameters["basic"]["MW"]
self.T_star = parameters["basic"]["T_star"]
self.rho_star = parameters["basic"]["rho_star"]
self.Tc = parameters["basic"]["Tc"]
self.rhoc = parameters["basic"]["rhoc"]
self.Pc = parameters["basic"]["Pc"]
self.Tt = parameters["basic"]["Tt"]
self.Pt = parameters["basic"]["Pt"]
self.rhot_l = parameters["basic"]["rhot_l"]
self.rhot_v = parameters["basic"]["rhot_v"]
self.P_min = parameters["basic"]["P_min"]
self.P_max = parameters["basic"]["P_max"]
self.rho_max = parameters["basic"]["rho_max"]
self.T_min = parameters["basic"]["T_min"]
self.T_max = parameters["basic"]["T_max"]
# If provided, set the reference state offset, otherwise set to (0, 0).
try:
self.reference_state_offset = parameters["eos"]["reference_state_offset"]
except KeyError:
self.reference_state_offset = [0.0, 0.0]
# Create the main Pyomo model to write equation of state expressions
self.model = self.make_model("delta", "tau")
# Optional models (Thermal Conductivity, Viscosity, and Surface Tension)
self.model_tcx = self.make_model("delta", "tau")
self.model_visc = self.make_model("delta", "tau")
# Since surface tension is only valid for the two phase region, surface tension
# is a function of temperature only
self.model_st = self.make_model("tau")
# Check if a predefined form of the ideal part of Helmholtz free energy is used
try:
phi_ideal_type = parameters["eos"]["phi_ideal_type"]
if phi_ideal_type:
self.add(
phi_ideal_types[phi_ideal_type](
model=self.model, parameters=parameters
)
)
except KeyError:
phi_ideal_type = 0
# Check if a predefined form of the residual part of Helmholtz free energy is used
try:
phi_residual_type = parameters["eos"]["phi_residual_type"]
if phi_residual_type:
self.add(
phi_residual_types[phi_residual_type](
model=self.model, parameters=parameters
)
)
except KeyError:
phi_residual_type = 0
# Check if predefined approximate liquid and vapor saturated density curves are used
aux_parameters = parameters.get("aux", None)
if aux_parameters is not None:
delta_l_sat_parameters = parameters["aux"].get("delta_l_sat_approx", None)
delta_v_sat_parameters = parameters["aux"].get("delta_v_sat_approx", None)
if delta_l_sat_parameters is not None:
etype = delta_l_sat_parameters.get("type", 0)
if etype:
self.add(
{
"delta_l_sat_approx": delta_sat_types[etype](
model=self.model,
name="delta_l_sat_approx",
parameters=parameters,
),
}
)
if delta_v_sat_parameters is not None:
etype = delta_v_sat_parameters.get("type", 0)
if etype:
self.add(
{
"delta_v_sat_approx": delta_sat_types[etype](
model=self.model,
name="delta_v_sat_approx",
parameters=parameters,
),
}
)
# Check if using predefined surface tension expression
try:
etype = parameters["transport"]["surface_tension"]["type"]
if etype:
self.add(
{
"surface_tension": surface_tension_types[etype](
model=self.model_st, parameters=parameters
),
}
)
except KeyError: # No surface tension to add
pass
[docs]
def calculate_pressure(self, rho, T):
"""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.
Args:
rho (float): density in kg/m3
T (float): temperature in K
Returns:
float: pressure in kPa
"""
self.model.delta = rho / self.rho_star
self.model.tau = self.T_star / T
return pyo.value(rho * self.R * T * (1 + self.model.delta * self.model.phir_d))
[docs]
def calculate_enthalpy(self, rho, T):
"""From the expressions provided, calculate enthalpy from density and
temperature. This can be used for testing.
Args:
rho (float): density in kg/m3
T (float): temperature in K
Returns:
float: enthalpy in kJ/kg
"""
self.model.delta = rho / self.rho_star
self.model.tau = self.T_star / T
return pyo.value(
self.R
* T
* (
1
+ self.model.tau * (self.model.phii_t + self.model.phir_t)
+ self.model.delta * self.model.phir_d
)
)
[docs]
def calculate_entropy(self, rho, T):
"""From the expressions provided, calculate entropy from density and
temperature. This can be used for testing.
Args:
rho (float): density in kg/m3
T (float): temperature in K
Returns:
float: entropy in kJ/kg/K
"""
self.model.delta = rho / self.rho_star
self.model.tau = self.T_star / T
return pyo.value(
self.R
* (
self.model.tau * (self.model.phii_t + self.model.phir_t)
- self.model.phii
- self.model.phir
)
)
[docs]
def make_model(self, *args):
"""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.
Args:
(str): Variables to add to the model
Returns:
ConcreteModel: Pyomo model with variables from args
"""
m = pyo.ConcreteModel()
for a in args:
setattr(m, a, pyo.Var())
m.R = pyo.Param(initialize=self.R)
m.MW = pyo.Param(initialize=self.MW)
m.T_star = pyo.Param(initialize=self.T_star)
m.rho_star = pyo.Param(initialize=self.rho_star)
m.Tc = pyo.Param(initialize=self.Tc)
m.rhoc = pyo.Param(initialize=self.rhoc)
m.Pc = pyo.Param(initialize=self.Pc)
m.Tt = pyo.Param(initialize=self.Tt)
m.Pt = pyo.Param(initialize=self.Pt)
m.rhot_l = pyo.Param(initialize=self.rhot_l)
m.rhot_v = pyo.Param(initialize=self.rhot_v)
m.P_min = pyo.Param(initialize=self.P_min)
m.P_max = pyo.Param(initialize=self.P_max)
m.rho_max = pyo.Param(initialize=self.rho_max)
m.T_min = pyo.Param(initialize=self.T_min)
m.T_max = pyo.Param(initialize=self.T_max)
return m
[docs]
def add(self, expressions):
"""This adds expressions to the to the object. These expressions are written
to NL files to be used by external functions.
Args:
expressions (dict): Dictionary where the key is an expression name and the
value is a Pyomo expression.
Returns:
None
"""
for name, expr in expressions.items():
if name in self.expressions: # check if in thermo model
m = self.model
elif name == "thermal_conductivity":
m = self.model_tcx
elif name == "surface_tension":
m = self.model_st
elif name == "viscosity":
m = self.model_visc
else:
raise RuntimeError(f"Unknown expression {name}")
if callable(expr):
setattr(m, name, pyo.Objective(rule=expr))
else:
setattr(m, name, pyo.Objective(expr=expr))
self.has_expression.append(name)
[docs]
def approx_sat_curves(self, trange):
"""_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.
Args:
trange (iterable): temperature points in K
Returns:
tuple: list of liquid densities and list of vapor densities
"""
_log.info(
"\n====================================================================="
)
_log.info(" Check approx sat delta curves")
_log.info(
"====================================================================="
)
_log.info(
f"{'T [K]':7s} {'T [C]':8s} {'~rho_l [kg/m3]':14s} {'~rho_v [kg/m3]':14s} {'~P [kPa]':9s}"
)
_log.info(
"---------------------------------------------------------------------"
)
rhol_list = []
rhov_list = []
for T in trange:
self.model.tau = self.model.T_star / T
delta_l = pyo.value(self.model.delta_l_sat_approx)
delta_v = pyo.value(self.model.delta_v_sat_approx)
rho_l = pyo.value(delta_l * self.model.rho_star)
rho_v = pyo.value(delta_v * self.model.rho_star)
rhol_list.append(rho_l)
rhov_list.append(rho_v)
self.model.delta = delta_v
P_v = pyo.value(
rho_v * self.R * T * (1 + self.model.delta * self.model.phir_d)
)
_log.info(
f"{T:7.3f}, {T - 273.15: 8.3f}, {rho_l:14.4f}, {rho_v:14.4f}, {P_v:9.4f}"
)
_log.info(
"=====================================================================\n"
)
return (rhol_list, rhov_list)
[docs]
def calculate_reference_offset(self, delta, tau, s0, h0):
"""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.
Args:
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:
(tuple): Offset for given reference state
"""
self.model.tau = tau
self.model.delta = delta
s1 = pyo.value(
self.model.tau * (self.model.phii_t + self.model.phir_t)
- self.model.phii
- self.model.phir
)
n1_off = s1 - s0
h1 = pyo.value(
1.0
/ self.model.tau
* (
1
+ self.model.tau * (self.model.phii_t + self.model.phir_t)
+ self.model.delta * self.model.phir_d
)
)
n2_off = h0 - h1
return (n1_off, n2_off)
[docs]
def write_model(self, model, model_name, expressions=None):
"""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.
Args:
model (ConcreteModel): a Pyomo model
model_name (str): the model name used in the NL file
Returns:
tuple: NL file, expression map, variable map for EOS
"""
nl_file, smap_id = model.write(f"{self.comp}_expressions_{model_name}.nl")
for v in model.component_data_objects(pyo.Var):
v.unfix()
smap = model.solutions.symbol_map[smap_id]
var_map = [1000] * 4
for s, c in smap.bySymbol.items():
if s.startswith("v"):
j = int(s[1:])
var_map[j] = self.variables[c.name]
if expressions is not None:
expr_map = [0] * len(expressions)
for s, c in smap.bySymbol.items():
if s.startswith("o"):
i = expressions[c.name]
j = int(s[1:])
expr_map[i] = j
return nl_file, expr_map, var_map
return nl_file, None, var_map
[docs]
def write(self, dry_run=False):
"""Write the parameter and expression files, and _log.info some diagnostics."""
if dry_run:
# TODO<jce> this isn't important to IDAES function, so I need to come back later
# and figure out how to test this.
return
_log.info(
"\n======================================================================="
)
_log.info(f" Writing expressions for {self.comp}")
_log.info(
"======================================================================="
)
pc = self.calculate_pressure(self.rhoc, self.Tc)
_log.info(f"Calculated Pc = {pc}, Pc given {self.Pc}")
_log.info(f"Tc = {self.Tc}, T^* = {self.T_star}")
_log.info(f"rhoc = {self.rhoc}, rho^* = {self.rho_star}")
self.Pc = pc
for name in self.expressions:
if name not in self.has_expression:
raise RuntimeError(f"Required expression {name} not provided.")
nl_file, expr_map, var_map = self.write_model(
self.model, "eos", self.expressions
)
param_dict = {
"nl_file": nl_file,
"expr_map": expr_map,
"var_map": var_map,
"param": {
"R": self.R,
"MW": self.MW,
"T_star": self.T_star,
"rho_star": self.rho_star,
"Tc": self.Tc,
"rhoc": self.rhoc,
"Pc": self.Pc,
"Tt": self.Tt,
"Pt": self.Pt,
"rhot_l": self.rhot_l,
"rhot_v": self.rhot_v,
"P_min": self.P_min,
"P_max": self.P_max,
"rho_max": self.rho_max,
"T_min": self.T_min,
"T_max": self.T_max,
"reference_state_offset": self.reference_state_offset,
},
}
# Add optional models
for name, short_name in self.optional_expressions.items():
if name in self.has_expression:
model = getattr(self, f"model_{short_name}")
nl_file, _, var_map = self.write_model(model, short_name)
param_dict[f"nl_file_{short_name}"] = nl_file
param_dict[f"var_map_{short_name}"] = var_map
param_dict[f"have_{short_name}"] = True
else:
_log.warning(f"Missing optional expression {name}")
param_dict[f"have_{short_name}"] = False
with open(f"{self.comp}_parameters.json", "w") as f:
json.dump(param_dict, f, indent=4)
trng = []
linspc = numpy.linspace(self.Tt, self.T_star, 15)
for i, t in enumerate(linspc):
if i != 0 and i != len(linspc) - 1:
t = round(t, 3)
trng.append(t)
self.approx_sat_curves(trng)
