#################################################################################
# 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.
#################################################################################
"""
General Cubic Equation of State property package with VLE calucations.
Cubic formulation and pure component property correlations from:
"The Properties of Gases and Liquids, 4th Edition", Reid, Prausnitz and Poling,
McGraw-Hill, 1987
Smooth Vapor-Liquid Equilibrium formulation from:
"A Smooth, Square Flash Formulation for Equation-Oriented Flowsheet
Optimization", Burgard et al., Proceedings of the 13 the International
Symposium on Process Systems Engineering – PSE 2018, July 1-5, 2018, San Diego
All results have been cross-referenced against other sources.
"""
# TODO: Missing docstrings
# pylint: disable=missing-function-docstring
# TODO: Look into protected access issues
# pylint: disable=protected-access
# Import Python libraries
import math
# Import Pyomo libraries
from pyomo.environ import (
Block,
check_optimal_termination,
Constraint,
exp,
Expression,
log,
NonNegativeReals,
sqrt,
Param,
PositiveReals,
value,
Var,
units as pyunits,
)
from pyomo.common.config import ConfigDict, ConfigValue, In
from pyomo.contrib.incidence_analysis import solve_strongly_connected_components
from pyomo.common.deprecation import deprecated
# Import IDAES cores
from idaes.core import (
declare_process_block_class,
MaterialFlowBasis,
PhysicalParameterBlock,
StateBlockData,
StateBlock,
MaterialBalanceType,
EnergyBalanceType,
LiquidPhase,
VaporPhase,
)
from idaes.core.util.initialization import (
solve_indexed_blocks,
fix_state_vars,
revert_state_vars,
)
from idaes.core.util.exceptions import (
BurntToast,
ConfigurationError,
InitializationError,
)
from idaes.core.util.model_statistics import (
degrees_of_freedom,
number_activated_equalities,
)
from idaes.core.util.math import safe_log
from idaes.core.solvers import get_solver
from idaes.core.util.constants import Constants as const
import idaes.logger as idaeslog
import idaes.core.util.scaling as iscale
from idaes.core.initialization.initializer_base import InitializerBase
# cubic_roots_available is used elsewhere
# pylint: disable=W0611
from idaes.models.properties.modular_properties.eos.ceos_common import (
EoS_param,
cubic_roots_available,
CubicThermoExpressions,
CubicType as CubicEoS,
)
# Set up logger
_log = idaeslog.getLogger(__name__)
[docs]
@deprecated(
msg="The standalone cubic property package has been deprecated in favor of the "
"cubic equation of state for the modular property framework. This class will be "
"removed in the May 2025 release.",
version="2.7.0",
)
@declare_process_block_class("CubicParameterBlock")
class CubicParameterData(PhysicalParameterBlock):
"""
General Property Parameter Block Class
"""
# Config block for the _IdealStateBlock
CONFIG = PhysicalParameterBlock.CONFIG()
CONFIG.declare(
"valid_phase",
ConfigValue(
default=("Vap", "Liq"),
domain=In(["Liq", "Vap", ("Vap", "Liq"), ("Liq", "Vap")]),
description="Flag indicating the valid phase",
doc="""Flag indicating the valid phase for a given set of
conditions, and thus corresponding constraints should be included,
**default** - ('Vap', 'Liq').
**Valid values:** {
**'Liq'** - Liquid only,
**'Vap'** - Vapor only,
**('Vap', 'Liq')** - Vapor-liquid equilibrium,
**('Liq', 'Vap')** - Vapor-liquid equilibrium,}""",
),
)
[docs]
def build(self):
"""
Callable method for Block construction.
"""
super(CubicParameterData, self).build()
self._state_block_class = CubicStateBlock
# Create Phase objects
if (
self.config.valid_phase == ("Liq", "Vap")
or self.config.valid_phase == ("Vap", "Liq")
or self.config.valid_phase == "Liq"
):
self.Liq = LiquidPhase()
if (
self.config.valid_phase == ("Liq", "Vap")
or self.config.valid_phase == ("Vap", "Liq")
or self.config.valid_phase == "Vap"
):
self.Vap = VaporPhase()
self.set_default_scaling("flow_mol", 1e-3)
self.set_default_scaling("mole_frac_comp", 10)
self.set_default_scaling("temperature", 1e-1)
self.set_default_scaling("temperature_dew", 1e-1)
self.set_default_scaling("temperature_bubble", 1e-1)
self.set_default_scaling("flow_mol_phase", 1e-2)
self.set_default_scaling("dens_mol_phase", 1)
self.set_default_scaling("dens_mass_phase", 1e-1)
self.set_default_scaling("pressure", 1e-6)
self.set_default_scaling("pressure_sat", 1e-6)
self.set_default_scaling("pressure_dew", 1e-6)
self.set_default_scaling("pressure_bubble", 1e-6)
self.set_default_scaling("mole_frac_phase_comp", 10)
self.set_default_scaling("enth_mol_phase", 1e-3, index="Liq")
self.set_default_scaling("enth_mol_phase", 1e-4, index="Vap")
self.set_default_scaling("enth_mol", 1e-3)
self.set_default_scaling("entr_mol_phase", 1e-2)
self.set_default_scaling("entr_mol", 1e-2)
self.set_default_scaling("mw", 100)
self.set_default_scaling("mw_phase", 100)
self.set_default_scaling("fug_phase_comp", 1)
self.set_default_scaling("fug_coeff_phase_comp", 1)
self.set_default_scaling("gibbs_mol_phase", 1e-3)
[docs]
@deprecated(
msg="The standalone cubic property package has been deprecated in favor of the "
"cubic equation of state for the modular property framework. This class will be "
"removed in the May 2025 release.",
version="2.7.0",
)
class CubicEoSInitializer(InitializerBase):
"""
Initializer for CubicEoS property packages.
This Initializer uses a sequential routine to initialize the
property package using the following steps:
1. Initialize bubble and dew point calculations (if present)
2. Estimate vapor-liquid equilibrium T_eq (if present)
3. Solve for phase-equilibrium conditions
4. Initialize all remaining properties
The Pyomo solve_strongly_connected_components method is used at each
step to converge the problem.
Note that for systems without vapor-liquid equilibrium the generic
BlockTriangularizationInitializer is probably sufficient for initializing
the property package.
"""
CONFIG = InitializerBase.CONFIG()
CONFIG.declare(
"solver",
ConfigValue(
default="ipopt_v2",
description="Solver to use for initialization",
),
)
CONFIG.declare(
"solver_options",
ConfigDict(
implicit=True,
description="Dict of options to pass to solver",
),
)
CONFIG.declare(
"solver_writer_config",
ConfigDict(
implicit=True,
description="Dict of writer_config arguments to pass to solver",
),
)
CONFIG.declare(
"calculate_variable_options",
ConfigDict(
implicit=True,
description="Dict of options to pass to 1x1 block solver",
doc="Dict of options to pass to calc_var_kwds argument in "
"solve_strongly_connected_components method. NOTE: models "
"involving ExternalFunctions must set "
"'diff_mode=differentiate.Modes.reverse_numeric'",
),
)
def __init__(self, **kwargs):
super().__init__(**kwargs)
self._solver = None
[docs]
def initialization_routine(
self,
model: Block,
):
"""
Sequential initialization routine for cubic EoS properties.
Args:
model: model to be initialized
Returns:
None
"""
# Setup loggers
init_log = idaeslog.getInitLogger(
model.name, self.config.output_level, tag="properties"
)
solve_log = idaeslog.getSolveLogger(
model.name, self.config.output_level, tag="properties"
)
# Create solver object
solver_obj = get_solver(
solver=self.config.solver,
solver_options=self.config.solver_options,
writer_config=self.config.solver_writer_config,
)
init_log.info("Starting initialization routine")
# Deactivate the constraints specific for outlet block i.e.
# when defined state is False
for k in model.values():
if k.config.defined_state is False:
k.sum_mole_frac_out.deactivate()
# ---------------------------------------------------------------------
# If present, initialize bubble and dew point calculations
# Antoine equation
def antoine_P(b, j, T):
return pyunits.convert_value(
value(
10
** (
b.params.antoine_coeff_A[j]
- b.params.antoine_coeff_B[j]
/ (T + b.params.antoine_coeff_C[j])
)
),
from_units=pyunits.bar,
to_units=pyunits.Pa,
)
for k in model.values():
# Bubble temperature initialization
if hasattr(k, "_mole_frac_tbub"):
Tbub0 = 0
for j in k.params.component_list:
Tbub0 += value(
k.mole_frac_comp[j]
* (
k.params.antoine_coeff_B[j]
/ (
k.params.antoine_coeff_A[j]
- math.log10(
value(
pyunits.convert(
k.pressure, to_units=pyunits.bar
)
)
)
)
- k.params.antoine_coeff_C[j]
)
)
err = 1
counter = 0
while err > 1e-2 and counter < 100:
f = value(
sum(
antoine_P(k, j, Tbub0) * k.mole_frac_comp[j]
for j in k.params.component_list
)
- k.pressure
)
df = value(
sum(
k.mole_frac_comp[j]
* k.params.antoine_coeff_B[j]
* math.log(10)
* antoine_P(k, j, Tbub0)
/ (Tbub0 + k.params.antoine_coeff_C[j]) ** 2
for j in k.params.component_list
)
)
if f / df > 20:
Tbub1 = Tbub0 - 20
elif f / df < -20:
Tbub1 = Tbub0 + 20
else:
Tbub1 = Tbub0 - f / df
err = abs(Tbub1 - Tbub0)
Tbub0 = Tbub1
counter += 1
k.temperature_bubble.value = Tbub0
for j in k.params.component_list:
k._mole_frac_tbub[j].value = value(
k.mole_frac_comp[j] * antoine_P(k, j, Tbub0) / k.pressure
)
# Dew temperature initialization
if hasattr(k, "_mole_frac_tdew"):
Tdew0 = 0
for j in k.params.component_list:
Tdew0 += value(
k.mole_frac_comp[j]
* (
k.params.antoine_coeff_B[j]
/ (
k.params.antoine_coeff_A[j]
- math.log10(
value(
pyunits.convert(
k.pressure, to_units=pyunits.bar
)
)
)
)
- k.params.antoine_coeff_C[j]
)
)
err = 1
counter = 0
while err > 1e-2 and counter < 100:
f = value(
k.pressure
* sum(
k.mole_frac_comp[j] / antoine_P(k, j, Tdew0)
for j in k.params.component_list
)
- 1
)
df = -value(
k.pressure
* math.log(10)
* sum(
k.mole_frac_comp[j]
* k.params.antoine_coeff_B[j]
/ (
(Tdew0 + k.params.antoine_coeff_C[j]) ** 2
* antoine_P(k, j, Tdew0)
)
for j in k.params.component_list
)
)
if f / df > 20:
Tdew1 = Tdew0 - 20
elif f / df < -20:
Tdew1 = Tdew0 + 20
else:
Tdew1 = Tdew0 - f / df
err = abs(Tdew1 - Tdew0)
Tdew0 = Tdew1
counter += 1
k.temperature_dew.value = Tdew0
for j in k.params.component_list:
k._mole_frac_tdew[j].value = value(
k.mole_frac_comp[j] * k.pressure / antoine_P(k, j, Tdew0)
)
# Bubble pressure initialization
if hasattr(k, "_mole_frac_pbub"):
k.pressure_bubble.value = value(
sum(
k.mole_frac_comp[j] * antoine_P(k, j, k.temperature)
for j in k.params.component_list
)
)
for j in k.params.component_list:
k._mole_frac_pbub[j].value = value(
k.mole_frac_comp[j]
* antoine_P(k, j, k.temperature)
/ k.pressure_bubble
)
k.pressure_bubble.display()
k._mole_frac_pbub.display()
# Dew pressure initialization
if hasattr(k, "_mole_frac_pdew"):
k.pressure_dew.value = value(
sum(
1 / (k.mole_frac_comp[j] / antoine_P(k, j, k.temperature))
for j in k.params.component_list
)
)
for j in k.params.component_list:
k._mole_frac_pdew[j].value = value(
k.mole_frac_comp[j]
* k.pressure_bubble
/ antoine_P(k, j, k.temperature)
)
# Solve bubble and dew point constraints
for c in k.component_objects(Constraint):
# Deactivate all property constraints
if c.local_name not in (
"eq_pressure_dew",
"eq_pressure_bubble",
"eq_temperature_dew",
"eq_temperature_bubble",
"_sum_mole_frac_tbub",
"_sum_mole_frac_tdew",
"_sum_mole_frac_pbub",
"_sum_mole_frac_pdew",
):
c.deactivate()
if number_activated_equalities(k) > 0:
with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
solve_strongly_connected_components(
k,
solver=solver_obj,
solve_kwds={"tee": slc.tee},
calc_var_kwds=self.config.calculate_variable_options,
)
init_log.info("Dew and bubble point initialization complete.")
# ---------------------------------------------------------------------
# If flash, initialize T1 and Teq
for k in model.values():
if (k.params.config.valid_phase == ("Liq", "Vap")) or (
k.params.config.valid_phase == ("Vap", "Liq")
):
k._t1.value = max(k.temperature.value, k.temperature_bubble.value)
k._teq.value = min(k._t1.value, k.temperature_dew.value)
init_log.info("Equilibrium temperature initialization complete.")
# ---------------------------------------------------------------------
# Initialize flow rates and compositions
for k in model.values():
if k.params.config.valid_phase == "Liq":
k.flow_mol_phase["Liq"].value = k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Liq", j].value = k.mole_frac_comp[j].value
elif k.params.config.valid_phase == "Vap":
k.flow_mol_phase["Vap"].value = k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Vap", j].value = k.mole_frac_comp[j].value
else:
if k.temperature.value > k.temperature_dew.value:
# Pure vapour
k.flow_mol_phase["Vap"].value = k.flow_mol.value
k.flow_mol_phase["Liq"].value = 1e-5 * k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Vap", j].value = k.mole_frac_comp[
j
].value
k.mole_frac_phase_comp["Liq", j].value = k._mole_frac_tdew[
j
].value
elif k.temperature.value < k.temperature_bubble.value:
# Pure liquid
k.flow_mol_phase["Vap"].value = 1e-5 * k.flow_mol.value
k.flow_mol_phase["Liq"].value = k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Vap", j].value = k._mole_frac_tbub[
j
].value
k.mole_frac_phase_comp["Liq", j].value = k.mole_frac_comp[
j
].value
else:
# Two-phase
# Estimate vapor fraction from distance from Tbub and Tdew
# Thanks to Rahul Gandhi for the method
vapRatio = value(
(k.temperature - k.temperature_bubble)
/ (k.temperature_dew - k.temperature_bubble)
)
k.flow_mol_phase["Vap"].value = value(vapRatio * k.flow_mol)
k.flow_mol_phase["Liq"].value = value((1 - vapRatio) * k.flow_mol)
# Initialize compositions using Rachford-Rice equation
for j in k.params.component_list:
kfact = value(antoine_P(k, j, k.temperature.value) / k.pressure)
k.mole_frac_phase_comp["Liq", j].value = value(
k.mole_frac_comp[j] / (1 + vapRatio * (kfact - 1))
)
k.mole_frac_phase_comp["Vap", j].value = value(
k.mole_frac_phase_comp["Liq", j] * kfact
)
# ---------------------------------------------------------------------
# Solve phase equilibrium constraints
for k in model.values():
for c in k.component_objects(Constraint):
# Activate equilibrium constraints
if c.local_name in (
"total_flow_balance",
"component_flow_balances",
"equilibrium_constraint",
"sum_mole_frac",
"_t1_constraint",
"_teq_constraint",
):
c.activate()
with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
solve_strongly_connected_components(
k,
solver=solver_obj,
solve_kwds={"tee": slc.tee},
calc_var_kwds=self.config.calculate_variable_options,
)
init_log.info("Phase equilibrium initialization complete")
# ---------------------------------------------------------------------
# Initialize other properties
for k in model.values():
for c in k.component_objects(Constraint):
# Activate all constraints except sum_mole_frac_out
if c.local_name not in ("sum_mole_frac_out"):
c.activate()
with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
solve_strongly_connected_components(
k,
solver=solver_obj,
solve_kwds={"tee": slc.tee},
calc_var_kwds=self.config.calculate_variable_options,
)
init_log.info("Property initialization routine finished.")
return None
class _CubicStateBlock(StateBlock):
"""
This Class contains methods which should be applied to Property Blocks as a
whole, rather than individual elements of indexed Property Blocks.
"""
# Set default initializer
default_initializer = CubicEoSInitializer
def fix_initialization_states(self):
"""
Fixes state variables for state blocks.
Returns:
None
"""
# Fix state variables
fix_state_vars(self)
# Also need to deactivate sum of mole fraction constraint
for k in self.values():
k.sum_mol_frac_out.deactivate()
def initialize(
blk,
state_args=None,
state_vars_fixed=False,
hold_state=False,
outlvl=idaeslog.NOTSET,
solver=None,
optarg=None,
):
"""
Initialization routine for property package.
Keyword Arguments:
state_args : Dictionary with initial guesses for the state vars
chosen. Note that if this method is triggered
through the control volume, and if initial guesses
were not provided at the unit model level, the
control volume passes the inlet values as initial
guess. Expected keys in state_args dict are:
* flow_mol
* mole_frac_comp (dict with components as keys)
* pressure
* temperature
outlvl : sets output level of initialization routine
optarg : solver options dictionary object (default=None, use
default solver options)
state_vars_fixed: Flag to denote if state vars have already been
fixed.
- True - states have already been fixed and
initialization does not need to worry
about fixing and unfixing variables.
- False - states have not been fixed. The state
block will deal with fixing/unfixing.
solver : str indicating which solver to use during
initialization (default = None, use default solver)
hold_state : flag indicating whether the initialization routine
should unfix any state variables fixed during
initialization (default=False).
- True - states variables are not unfixed, and
a dict of returned containing flags for
which states were fixed during
initialization.
- False - state variables are unfixed after
initialization by calling the
release_state method
Returns:
If hold_states is True, returns a dict containing flags for
which states were fixed during initialization.
"""
init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties")
solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties")
init_log.info("Starting initialization")
# Deactivate the constraints specific for outlet block i.e.
# when defined state is False
for k in blk.values():
if k.config.defined_state is False:
k.sum_mole_frac_out.deactivate()
# Fix state variables if not already fixed
if state_vars_fixed is False:
flags = fix_state_vars(blk, state_args)
else:
# Check when the state vars are fixed already result in dof 0
for k in blk.values():
if degrees_of_freedom(k) != 0:
# PYLINT-TODO
# pylint: disable-next=broad-exception-raised
raise Exception(
"State vars fixed but degrees of freedom "
"for state block is not zero during "
"initialization."
)
# Create solver
opt = get_solver(solver, optarg)
# ---------------------------------------------------------------------
# If present, initialize bubble and dew point calculations
# Antoine equation
def antoine_P(b, j, T):
return pyunits.convert_value(
value(
10
** (
b.params.antoine_coeff_A[j]
- b.params.antoine_coeff_B[j]
/ (T + b.params.antoine_coeff_C[j])
)
),
from_units=pyunits.bar,
to_units=pyunits.Pa,
)
# Bubble temperature initialization
for k in blk.values():
if hasattr(k, "_mole_frac_tbub"):
Tbub0 = 0
for j in k.params.component_list:
Tbub0 += value(
k.mole_frac_comp[j]
* (
k.params.antoine_coeff_B[j]
/ (
k.params.antoine_coeff_A[j]
- math.log10(
value(
pyunits.convert(
k.pressure, to_units=pyunits.bar
)
)
)
)
- k.params.antoine_coeff_C[j]
)
)
err = 1
counter = 0
while err > 1e-2 and counter < 100:
f = value(
sum(
antoine_P(k, j, Tbub0) * k.mole_frac_comp[j]
for j in k.params.component_list
)
- k.pressure
)
df = value(
sum(
k.mole_frac_comp[j]
* k.params.antoine_coeff_B[j]
* math.log(10)
* antoine_P(k, j, Tbub0)
/ (Tbub0 + k.params.antoine_coeff_C[j]) ** 2
for j in k.params.component_list
)
)
if f / df > 20:
Tbub1 = Tbub0 - 20
elif f / df < -20:
Tbub1 = Tbub0 + 20
else:
Tbub1 = Tbub0 - f / df
err = abs(Tbub1 - Tbub0)
Tbub0 = Tbub1
counter += 1
k.temperature_bubble.value = Tbub0
for j in k.params.component_list:
k._mole_frac_tbub[j].value = value(
k.mole_frac_comp[j] * antoine_P(k, j, Tbub0) / k.pressure
)
# Dew temperature initialization
for k in blk.values():
if hasattr(k, "_mole_frac_tdew"):
Tdew0 = 0
for j in k.params.component_list:
Tdew0 += value(
k.mole_frac_comp[j]
* (
k.params.antoine_coeff_B[j]
/ (
k.params.antoine_coeff_A[j]
- math.log10(
value(
pyunits.convert(
k.pressure, to_units=pyunits.bar
)
)
)
)
- k.params.antoine_coeff_C[j]
)
)
err = 1
counter = 0
while err > 1e-2 and counter < 100:
f = value(
k.pressure
* sum(
k.mole_frac_comp[j] / antoine_P(k, j, Tdew0)
for j in k.params.component_list
)
- 1
)
df = -value(
k.pressure
* math.log(10)
* sum(
k.mole_frac_comp[j]
* k.params.antoine_coeff_B[j]
/ (
(Tdew0 + k.params.antoine_coeff_C[j]) ** 2
* antoine_P(k, j, Tdew0)
)
for j in k.params.component_list
)
)
if f / df > 20:
Tdew1 = Tdew0 - 20
elif f / df < -20:
Tdew1 = Tdew0 + 20
else:
Tdew1 = Tdew0 - f / df
err = abs(Tdew1 - Tdew0)
Tdew0 = Tdew1
counter += 1
k.temperature_dew.value = Tdew0
for j in k.params.component_list:
k._mole_frac_tdew[j].value = value(
k.mole_frac_comp[j] * k.pressure / antoine_P(k, j, Tdew0)
)
# Bubble pressure initialization
for k in blk.values():
if hasattr(k, "_mole_frac_pbub"):
k.pressure_bubble.value = value(
sum(
k.mole_frac_comp[j] * antoine_P(k, j, k.temperature)
for j in k.params.component_list
)
)
for j in k.params.component_list:
k._mole_frac_pbub[j].value = value(
k.mole_frac_comp[j]
* antoine_P(k, j, k.temperature)
/ k.pressure_bubble
)
k.pressure_bubble.display()
k._mole_frac_pbub.display()
# Dew pressure initialization
for k in blk.values():
if hasattr(k, "_mole_frac_pdew"):
k.pressure_dew.value = value(
sum(
1 / (k.mole_frac_comp[j] / antoine_P(k, j, k.temperature))
for j in k.params.component_list
)
)
for j in k.params.component_list:
k._mole_frac_pdew[j].value = value(
k.mole_frac_comp[j]
* k.pressure_bubble
/ antoine_P(k, j, k.temperature)
)
# Solve bubble and dew point constraints
cons_count = 0
for k in blk.values():
for c in k.component_objects(Constraint):
# Deactivate all property constraints
if c.local_name not in (
"eq_pressure_dew",
"eq_pressure_bubble",
"eq_temperature_dew",
"eq_temperature_bubble",
"_sum_mole_frac_tbub",
"_sum_mole_frac_tdew",
"_sum_mole_frac_pbub",
"_sum_mole_frac_pdew",
):
c.deactivate()
cons_count += number_activated_equalities(k)
if cons_count > 0:
with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
else:
res = ""
init_log.info(
"Dew and bubble point init complete {}.".format(idaeslog.condition(res))
)
# ---------------------------------------------------------------------
# If flash, initialize T1 and Teq
for k in blk.values():
if (k.params.config.valid_phase == ("Liq", "Vap")) or (
k.params.config.valid_phase == ("Vap", "Liq")
):
k._t1.value = max(k.temperature.value, k.temperature_bubble.value)
k._teq.value = min(k._t1.value, k.temperature_dew.value)
init_log.info("Equilibrium temperature init complete.")
# ---------------------------------------------------------------------
# Initialize flow rates and compositions
# TODO : This will need to be generalised more when we move to a
# modular implementation
for k in blk.values():
if k.params.config.valid_phase == "Liq":
k.flow_mol_phase["Liq"].value = k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Liq", j].value = k.mole_frac_comp[j].value
elif k.params.config.valid_phase == "Vap":
k.flow_mol_phase["Vap"].value = k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Vap", j].value = k.mole_frac_comp[j].value
else:
if k.temperature.value > k.temperature_dew.value:
# Pure vapour
k.flow_mol_phase["Vap"].value = k.flow_mol.value
k.flow_mol_phase["Liq"].value = 1e-5 * k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Vap", j].value = k.mole_frac_comp[
j
].value
k.mole_frac_phase_comp["Liq", j].value = k._mole_frac_tdew[
j
].value
elif k.temperature.value < k.temperature_bubble.value:
# Pure liquid
k.flow_mol_phase["Vap"].value = 1e-5 * k.flow_mol.value
k.flow_mol_phase["Liq"].value = k.flow_mol.value
for j in k.params.component_list:
k.mole_frac_phase_comp["Vap", j].value = k._mole_frac_tbub[
j
].value
k.mole_frac_phase_comp["Liq", j].value = k.mole_frac_comp[
j
].value
else:
# Two-phase
# Estimate vapor fraction from distance from Tbub and Tdew
# Thanks to Rahul Gandhi for the method
vapRatio = value(
(k.temperature - k.temperature_bubble)
/ (k.temperature_dew - k.temperature_bubble)
)
k.flow_mol_phase["Vap"].value = value(vapRatio * k.flow_mol)
k.flow_mol_phase["Liq"].value = value((1 - vapRatio) * k.flow_mol)
# Initialize compositions using Rachford-Rice equation
for j in k.params.component_list:
kfact = value(antoine_P(k, j, k.temperature.value) / k.pressure)
k.mole_frac_phase_comp["Liq", j].value = value(
k.mole_frac_comp[j] / (1 + vapRatio * (kfact - 1))
)
k.mole_frac_phase_comp["Vap", j].value = value(
k.mole_frac_phase_comp["Liq", j] * kfact
)
# ---------------------------------------------------------------------
# Solve phase equilibrium constraints
for k in blk.values():
for c in k.component_objects(Constraint):
# Activate equilibrium constraints
if c.local_name in (
"total_flow_balance",
"component_flow_balances",
"equilibrium_constraint",
"sum_mole_frac",
"_t1_constraint",
"_teq_constraint",
):
c.activate()
with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
results = solve_indexed_blocks(opt, [blk], tee=slc.tee)
init_log.info("Phase equilibrium init: {}.".format(idaeslog.condition(results)))
# ---------------------------------------------------------------------
# Initialize other properties
for k in blk.values():
for c in k.component_objects(Constraint):
# Activate all constraints except sum_mole_frac_out
if c.local_name not in ("sum_mole_frac_out"):
c.activate()
with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
results = solve_indexed_blocks(opt, [blk], tee=slc.tee)
init_log.info("Property init: {}.".format(idaeslog.condition(results)))
# ---------------------------------------------------------------------
if not check_optimal_termination(results):
raise InitializationError(
f"{blk.name} failed to initialize successfully. Please check "
f"the output logs for more information."
)
if state_vars_fixed is False:
if hold_state is True:
return flags
else:
blk.release_state(flags, outlvl=outlvl)
init_log.info("Initialization complete.")
def release_state(blk, flags, outlvl=idaeslog.NOTSET):
"""
Method to release state variables fixed during initialization.
Keyword Arguments:
flags : dict containing information of which state variables
were fixed during initialization, and should now be
unfixed. This dict is returned by initialize if
hold_state=True.
outlvl : sets output level of of logging
"""
for k in blk.values():
if not k.config.defined_state:
k.sum_mole_frac_out.activate()
if flags is None:
return
# Unfix state variables
revert_state_vars(blk, flags)
init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties")
init_log.info_high("States released.")
[docs]
@declare_process_block_class("CubicStateBlock", block_class=_CubicStateBlock)
class CubicStateBlockData(StateBlockData):
"""An general property package for cubic equations of state with VLE."""
[docs]
def build(self):
"""Callable method for Block construction."""
super(CubicStateBlockData, self).build()
self.expr_write = CubicThermoExpressions(self)
# Add state variables
self.flow_mol = Var(
initialize=1.0,
domain=NonNegativeReals,
doc="Component molar flowrate [mol/s]",
units=pyunits.mol / pyunits.s,
)
self.mole_frac_comp = Var(
self.params.component_list,
bounds=(0, None),
initialize=1 / len(self.params.component_list),
doc="Mixture mole fractions [-]",
)
self.pressure = Var(
initialize=101325,
domain=NonNegativeReals,
doc="State pressure [Pa]",
units=pyunits.Pa,
)
self.temperature = Var(
initialize=298.15,
domain=NonNegativeReals,
doc="State temperature [K]",
units=pyunits.K,
)
# Add supporting variables
self.flow_mol_phase = Var(
self.params.phase_list,
initialize=0.5,
domain=NonNegativeReals,
doc="Phase molar flow rates [mol/s]",
units=pyunits.mol / pyunits.s,
)
self.mole_frac_phase_comp = Var(
self.params.phase_list,
self.params.component_list,
initialize=1 / len(self.params.component_list),
bounds=(0, None),
doc="Phase mole fractions [-]",
)
if self.params.config.valid_phase == "Liq":
self._make_liq_phase_eq()
elif self.params.config.valid_phase == "Vap":
self._make_vap_phase_eq()
elif (self.params.config.valid_phase == ("Liq", "Vap")) or (
self.params.config.valid_phase == ("Vap", "Liq")
):
self._make_flash_eq()
else:
raise BurntToast(
"{} found unexpected value for valid_phases. "
"Please contact the "
"IDAES developers with this bug.".format(self.name)
)
def _make_liq_phase_eq(self):
# Add equilibrium temperature - in this case the state temperature
self._teq = Expression(expr=self.temperature)
# Add supporting equations for Cubic EoS
self.common_cubic()
def rule_total_mass_balance(b):
return b.flow_mol_phase["Liq"] == b.flow_mol
self.total_flow_balance = Constraint(rule=rule_total_mass_balance)
def rule_comp_mass_balance(b, i):
return b.mole_frac_comp[i] == b.mole_frac_phase_comp["Liq", i]
self.component_flow_balances = Constraint(
self.params.component_list, rule=rule_comp_mass_balance
)
if self.config.defined_state is False:
# applied at outlet only
self.sum_mole_frac_out = Constraint(
expr=1
== sum(self.mole_frac_comp[i] for i in self.params.component_list)
)
def _make_vap_phase_eq(self):
# Add equilibrium temperature - in this case the state temperature
self._teq = Expression(expr=self.temperature)
# Add supporting equations for Cubic EoS
self.common_cubic()
def rule_total_mass_balance(b):
return b.flow_mol_phase["Vap"] == b.flow_mol
self.total_flow_balance = Constraint(rule=rule_total_mass_balance)
def rule_comp_mass_balance(b, i):
return b.mole_frac_comp[i] == b.mole_frac_phase_comp["Vap", i]
self.component_flow_balances = Constraint(
self.params.component_list, rule=rule_comp_mass_balance
)
if self.config.defined_state is False:
# applied at outlet only
self.sum_mole_frac_out = Constraint(
expr=1
== sum(self.mole_frac_comp[i] for i in self.params.component_list)
)
def _make_flash_eq(self):
"""
Implementation of smooth VLE formulation.
See module header for reference.
"""
def rule_total_mass_balance(b):
return b.flow_mol_phase["Liq"] + b.flow_mol_phase["Vap"] == b.flow_mol
self.total_flow_balance = Constraint(rule=rule_total_mass_balance)
def rule_comp_mass_balance(b, i):
return (
b.flow_mol * b.mole_frac_comp[i]
== b.flow_mol_phase["Liq"] * b.mole_frac_phase_comp["Liq", i]
+ b.flow_mol_phase["Vap"] * b.mole_frac_phase_comp["Vap", i]
)
self.component_flow_balances = Constraint(
self.params.component_list, rule=rule_comp_mass_balance
)
def rule_mole_frac(b):
return (
sum(b.mole_frac_phase_comp["Liq", i] for i in b.params.component_list)
- sum(b.mole_frac_phase_comp["Vap", i] for i in b.params.component_list)
== 0
)
self.sum_mole_frac = Constraint(rule=rule_mole_frac)
if self.config.defined_state is False:
# applied at outlet only
self.sum_mole_frac_out = Constraint(
expr=1
== sum(self.mole_frac_comp[i] for i in self.params.component_list)
)
# Definition of equilibrium temperature for smooth VLE
self._teq = Var(
initialize=self.temperature.value,
doc="Temperature for calculating phase equilibrium",
units=pyunits.K,
)
self._t1 = Var(
initialize=self.temperature.value,
doc="Intermediate temperature for calculating Teq",
units=pyunits.K,
)
self.eps_1 = Param(
default=0.01,
mutable=True,
doc="Smoothing parameter for Teq",
units=pyunits.K,
)
self.eps_2 = Param(
default=0.0005,
mutable=True,
doc="Smoothing parameter for Teq",
units=pyunits.K,
)
# Add supporting equations for Cubic EoS
self.common_cubic()
# PSE paper Eqn 13
def rule_t1(b):
return b._t1 == 0.5 * (
b.temperature
+ b.temperature_bubble
+ sqrt((b.temperature - b.temperature_bubble) ** 2 + b.eps_1**2)
)
self._t1_constraint = Constraint(rule=rule_t1)
# PSE paper Eqn 14
# TODO : Add option for supercritical extension
def rule_teq(b):
return b._teq == 0.5 * (
b._t1
+ b.temperature_dew
- sqrt((b._t1 - b.temperature_dew) ** 2 + b.eps_2**2)
)
self._teq_constraint = Constraint(rule=rule_teq)
def rule_tr_eq(b, i):
return b._teq / b.params.temperature_crit[i]
self._tr_eq = Expression(
self.params.component_list,
rule=rule_tr_eq,
doc="Component reduced temperatures [-]",
)
def rule_equilibrium(b, i):
return b._log_equilibrium_cubic("Vap", i) == b._log_equilibrium_cubic(
"Liq", i
)
self.equilibrium_constraint = Constraint(
self.params.component_list, rule=rule_equilibrium
)
# -----------------------------------------------------------------------------
# Property Methods
def _dens_mol_phase(self):
self.dens_mol_phase = Var(
self.params.phase_list,
doc="Molar density [mol/m^3]",
units=pyunits.mol / pyunits.m**3,
)
def rule_dens_mol_phase(b, p):
if p == "Vap":
return b._dens_mol_vap()
else:
return b._dens_mol_liq()
self.eq_dens_mol_phase = Constraint(
self.params.phase_list, rule=rule_dens_mol_phase
)
def _dens_mass_phase(self):
self.dens_mass_phase = Var(
self.params.phase_list,
doc="Mass density [kg/m^3]",
units=pyunits.kg / pyunits.m**3,
)
def rule_dens_mass_phase(b, p):
if p == "Vap":
return b._dens_mass_vap()
else:
return b._dens_mass_liq()
self.eq_dens_mass_phase = Constraint(
self.params.phase_list, rule=rule_dens_mass_phase
)
def _energy_internal_mol_phase(self):
self.energy_internal_mol_phase = Var(
self.params.phase_list,
doc="Phase molar specific internal energies [J/mol]",
units=pyunits.J / pyunits.mol,
)
def rule_energy_internal_mol_phase(b, p):
if p == "Vap":
return b.energy_internal_mol_phase[p] == b._energy_internal_mol_vap()
else:
return b.energy_internal_mol_phase[p] == b._energy_internal_mol_liq()
self.eq_energy_internal_mol_phase = Constraint(
self.params.phase_list, rule=rule_energy_internal_mol_phase
)
def _energy_internal_mol(self):
self.energy_internal_mol = Var(
doc="Mixture molar specific internal energy [J/mol]",
units=pyunits.J / pyunits.mol,
)
def rule_energy_internal_mol(b):
return b.energy_internal_mol * b.flow_mol == sum(
b.flow_mol_phase[p] * b.energy_internal_mol_phase[p]
for p in b.params.phase_list
)
self.eq_energy_internal_mol = Constraint(rule=rule_energy_internal_mol)
def _enth_mol_phase(self):
self.enth_mol_phase = Var(
self.params.phase_list,
doc="Phase molar specific enthalpies [J/mol]",
units=pyunits.J / pyunits.mol,
)
def rule_enth_mol_phase(b, p):
if p == "Vap":
return b.enth_mol_phase[p] == b._enth_mol_vap()
else:
return b.enth_mol_phase[p] == b._enth_mol_liq()
self.eq_enth_mol_phase = Constraint(
self.params.phase_list, rule=rule_enth_mol_phase
)
def _enth_mol(self):
self.enth_mol = Var(
doc="Mixture molar specific enthalpies [J/mol]",
units=pyunits.J / pyunits.mol,
)
def rule_enth_mol(b):
return b.enth_mol * b.flow_mol == sum(
b.flow_mol_phase[p] * b.enth_mol_phase[p] for p in b.params.phase_list
)
self.eq_enth_mol = Constraint(rule=rule_enth_mol)
def _entr_mol(self):
self.entr_mol = Var(
doc="Mixture molar specific entropies [J/mol.K]",
units=pyunits.J / pyunits.mol / pyunits.K,
)
def rule_entr_mol(b):
return b.entr_mol * b.flow_mol == sum(
b.flow_mol_phase[p] * b.entr_mol_phase[p] for p in b.params.phase_list
)
self.eq_entr_mol = Constraint(rule=rule_entr_mol)
def _entr_mol_phase(self):
self.entr_mol_phase = Var(
self.params.phase_list,
doc="Phase molar specific entropies [J/mol.K]",
units=pyunits.J / pyunits.mol / pyunits.K,
)
def rule_entr_mol_phase(b, p):
if p == "Vap":
return b.entr_mol_phase[p] == b._entr_mol_vap()
else:
return b.entr_mol_phase[p] == b._entr_mol_liq()
self.eq_entr_mol_phase = Constraint(
self.params.phase_list, rule=rule_entr_mol_phase
)
def _fug_phase(self):
def rule_fug_phase(b, p, j):
if p == "Vap":
return b._fug_vap(j)
else:
return b._fug_liq(j)
self.fug_phase_comp = Expression(
self.params.phase_list, self.params.component_list, rule=rule_fug_phase
)
def _fug_coeff_phase(self):
def rule_fug_coeff_phase(b, p, j):
if p == "Vap":
return b._fug_coeff_vap(j)
else:
return b._fug_coeff_liq(j)
self.fug_coeff_phase_comp = Expression(
self.params.phase_list,
self.params.component_list,
rule=rule_fug_coeff_phase,
)
def _gibbs_mol_phase(self):
self.gibbs_mol_phase = Var(
self.params.phase_list,
doc="Phase molar specific Gibbs energy [J/mol]",
units=pyunits.J / pyunits.mol,
)
def rule_gibbs_mol_phase(b, p):
return b.gibbs_mol_phase[p] == (
b.enth_mol_phase[p] - b.temperature * b.entr_mol_phase[p]
)
self.eq_gibbs_mol_phase = Constraint(
self.params.phase_list, rule=rule_gibbs_mol_phase
)
def _mw(self):
def rule_mw(b):
return sum(b.mw_phase[p] for p in b.params.phase_list)
self.mw = Expression(rule=rule_mw)
def _mw_phase(self):
def rule_mw_phase(b, p):
return sum(
b.mole_frac_phase_comp[p, j] * b.params.mw_comp[j]
for j in b.params.component_list
)
self.mw_phase = Expression(self.params.phase_list, rule=rule_mw_phase)
# -----------------------------------------------------------------------------
# General Methods
[docs]
def get_material_flow_terms(self, p, j):
"""Create material flow terms for control volume."""
if not self.is_property_constructed("material_flow_terms"):
try:
def rule_material_flow_terms(b, p, j):
return self.flow_mol_phase[p] * self.mole_frac_phase_comp[p, j]
self.material_flow_terms = Expression(
self.params.phase_list,
self.params.component_list,
rule=rule_material_flow_terms,
)
except AttributeError:
self.del_component(self.material_flow_terms)
if j in self.params.component_list:
return self.material_flow_terms[p, j]
else:
return 0
[docs]
def get_enthalpy_flow_terms(self, p):
"""Create enthalpy flow terms."""
if not self.is_property_constructed("enthalpy_flow_terms"):
try:
def rule_enthalpy_flow_terms(b, p):
return self.flow_mol_phase[p] * self.enth_mol_phase[p]
self.enthalpy_flow_terms = Expression(
self.params.phase_list, rule=rule_enthalpy_flow_terms
)
except AttributeError:
self.del_component(self.enthalpy_flow_terms)
return self.enthalpy_flow_terms[p]
[docs]
def get_material_density_terms(self, p, j):
"""Create material density terms."""
if not self.is_property_constructed("material_density_terms"):
try:
def rule_material_density_terms(b, p, j):
return self.dens_mol_phase[p] * self.mole_frac_phase_comp[p, j]
self.material_density_terms = Expression(
self.params.phase_list,
self.params.component_list,
rule=rule_material_density_terms,
)
except AttributeError:
self.del_component(self.material_density_terms)
if j in self.params.component_list:
return self.material_density_terms[p, j]
else:
return 0
[docs]
def get_energy_density_terms(self, p):
"""Create energy density terms."""
if not self.is_property_constructed("energy_density_terms"):
try:
def rule_energy_density_terms(b, p):
return self.dens_mol_phase[p] * self.energy_internal_mol_phase[p]
self.energy_density_terms = Expression(
self.params.phase_list, rule=rule_energy_density_terms
)
except AttributeError:
self.del_component(self.energy_density_terms)
return self.energy_density_terms[p]
def default_material_balance_type(self):
return MaterialBalanceType.componentTotal
def default_energy_balance_type(self):
return EnergyBalanceType.enthalpyTotal
[docs]
def get_material_flow_basis(b):
return MaterialFlowBasis.molar
[docs]
def define_state_vars(self):
"""Define state vars."""
return {
"flow_mol": self.flow_mol,
"mole_frac_comp": self.mole_frac_comp,
"temperature": self.temperature,
"pressure": self.pressure,
}
[docs]
def define_display_vars(b):
return {
"Molar Flowrate": b.flow_mol,
"Mole Fractions": b.mole_frac_comp,
"Temperature": b.temperature,
"Pressure": b.pressure,
}
[docs]
def model_check(blk):
"""Model checks for property block."""
# Check temperature bounds
if value(blk.temperature) < blk.temperature.lb:
_log.error("{} Temperature set below lower bound.".format(blk.name))
if value(blk.temperature) > blk.temperature.ub:
_log.error("{} Temperature set above upper bound.".format(blk.name))
# Check pressure bounds
if value(blk.pressure) < blk.pressure.lb:
_log.error("{} Pressure set below lower bound.".format(blk.name))
if value(blk.pressure) > blk.pressure.ub:
_log.error("{} Pressure set above upper bound.".format(blk.name))
# -----------------------------------------------------------------------------
# Bubble and Dew Points
def _temperature_bubble(self):
self.temperature_bubble = Var(
doc="Bubble point temperature (K)", units=pyunits.K
)
self._mole_frac_tbub = Var(
self.params.component_list,
initialize=1 / len(self.params.component_list),
bounds=(0, None),
doc="Vapor mole fractions at bubble point",
)
self._sum_mole_frac_tbub = Constraint(
expr=1e3
== 1e3 * sum(self._mole_frac_tbub[j] for j in self.params.component_list)
)
def rule_bubble_temp(b, j):
return log(b.mole_frac_comp[j]) + log(b.bubble_temp_liq(j)) == log(
b._mole_frac_tbub[j]
) + log(b.bubble_temp_vap(j))
self.eq_temperature_bubble = Constraint(
self.params.component_list, rule=rule_bubble_temp
)
def _temperature_dew(self):
self.temperature_dew = Var(doc="Dew point temperature (K)", units=pyunits.K)
self._mole_frac_tdew = Var(
self.params.component_list,
initialize=1 / len(self.params.component_list),
bounds=(0, None),
doc="Liquid mole fractions at dew point",
)
self._sum_mole_frac_tdew = Constraint(
expr=1e3
== 1e3 * sum(self._mole_frac_tdew[j] for j in self.params.component_list)
)
def rule_dew_temp(b, j):
return log(b._mole_frac_tdew[j]) + log(b.dew_temp_liq(j)) == log(
b.mole_frac_comp[j]
) + log(b.dew_temp_vap(j))
self.eq_temperature_dew = Constraint(
self.params.component_list, rule=rule_dew_temp
)
def _pressure_bubble(self):
self.pressure_bubble = Var(
domain=PositiveReals, doc="Bubble point pressure (Pa)", units=pyunits.Pa
)
self._mole_frac_pbub = Var(
self.params.component_list,
initialize=1 / len(self.params.component_list),
bounds=(0, None),
doc="Vapor mole fractions at bubble point",
)
self._sum_mole_frac_pbub = Constraint(
expr=1e3
== 1e3 * sum(self._mole_frac_pbub[j] for j in self.params.component_list)
)
def rule_bubble_pres(b, j):
return log(b.mole_frac_comp[j]) + log(b.bubble_pres_liq(j)) == log(
b._mole_frac_pbub[j]
) + log(b.bubble_pres_vap(j))
self.eq_pressure_bubble = Constraint(
self.params.component_list, rule=rule_bubble_pres
)
def _pressure_dew(self):
self.pressure_dew = Var(
domain=PositiveReals, doc="Dew point pressure (Pa)", units=pyunits.Pa
)
self._mole_frac_pdew = Var(
self.params.component_list,
initialize=1 / len(self.params.component_list),
bounds=(0, None),
doc="Liquid mole fractions at dew point",
)
self._sum_mole_frac_pdew = Constraint(
expr=1 == sum(self._mole_frac_pdew[j] for j in self.params.component_list)
)
def rule_dew_press(b, j):
return log(b._mole_frac_pdew[j]) + log(b.dew_press_liq(j)) == log(
b.mole_frac_comp[j]
) + log(b.dew_press_vap(j))
self.eq_pressure_dew = Constraint(
self.params.component_list, rule=rule_dew_press
)
# -----------------------------------------------------------------------------
# Liquid phase properties
def _vol_mol_liq(b):
return b._vol_mol_cubic("Liq")
def _dens_mol_liq(b):
return b._dens_mol_cubic("Liq")
def _dens_mass_liq(b):
return b._dens_mass_cubic("Liq")
def _fug_liq(b, j):
return b._fug_cubic("Liq", j)
def _fug_coeff_liq(b, j):
return b._fug_coeff_cubic("Liq", j)
def _energy_internal_mol_liq(b):
return b._energy_internal_mol_cubic("Liq")
def _enth_mol_liq(b):
return b._enth_mol_cubic("Liq")
def _enth_mol_liq_ig(b):
return b._enth_mol_ig("Liq")
def _entr_mol_liq(b):
return b._entr_mol_cubic("Liq")
def _entr_mol_liq_ig(b):
return b._entr_mol_ig("Liq")
def bubble_temp_liq(b, j):
def a(k):
return (
b.omegaA
* (
(const.gas_constant * b.params.temperature_crit[k]) ** 2
/ b.params.pressure_crit[k]
)
* (
(
1
+ b.fw[k]
* (
1
- sqrt(b.temperature_bubble / b.params.temperature_crit[k])
)
)
** 2
)
)
am = sum(
sum(
b.mole_frac_comp[i]
* b.mole_frac_comp[j]
* sqrt(a(i) * a(j))
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b.mole_frac_comp[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure / (const.gas_constant * b.temperature_bubble) ** 2
B = bm * b.pressure / (const.gas_constant * b.temperature_bubble)
delta = (
2
* sqrt(a(j))
/ am
* sum(
b.mole_frac_comp[i] * sqrt(a(i)) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_liq(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
def dew_temp_liq(b, j):
def a(k):
return (
b.omegaA
* (
(const.gas_constant * b.params.temperature_crit[k]) ** 2
/ b.params.pressure_crit[k]
)
* (
(
1
+ b.fw[k]
* (1 - sqrt(b.temperature_dew / b.params.temperature_crit[k]))
)
** 2
)
)
am = sum(
sum(
b._mole_frac_tdew[i]
* b._mole_frac_tdew[j]
* sqrt(a(i) * a(j))
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b._mole_frac_tdew[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure / (const.gas_constant * b.temperature_dew) ** 2
B = bm * b.pressure / (const.gas_constant * b.temperature_dew)
delta = (
2
* sqrt(a(j))
/ am
* sum(
b._mole_frac_tdew[i] * sqrt(a(i)) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_liq(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
def bubble_pres_liq(b, j):
am = sum(
sum(
b.mole_frac_comp[i]
* b.mole_frac_comp[j]
* sqrt(b.a[i] * b.a[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b.mole_frac_comp[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure_bubble / (const.gas_constant * b.temperature) ** 2
B = bm * b.pressure_bubble / (const.gas_constant * b.temperature)
delta = (
2
* sqrt(b.a[j])
/ am
* sum(
b.mole_frac_comp[i] * sqrt(b.a[i]) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_liq(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
def dew_press_liq(b, j):
am = sum(
sum(
b._mole_frac_pdew[i]
* b._mole_frac_pdew[j]
* sqrt(b.a[i] * b.a[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b._mole_frac_pdew[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure_dew / (const.gas_constant * b.temperature) ** 2
B = bm * b.pressure_dew / (const.gas_constant * b.temperature)
delta = (
2
* sqrt(b.a[j])
/ am
* sum(
b._mole_frac_pdew[i] * sqrt(b.a[i]) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_liq(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
# -----------------------------------------------------------------------------
# Vapour phase properties
def _vol_mol_vap(b):
return b._vol_mol_cubic("Vap")
def _dens_mol_vap(b):
return b._dens_mol_cubic("Vap")
def _dens_mass_vap(b):
return b._dens_mass_cubic("Vap")
def _fug_vap(b, j):
return b._fug_cubic("Vap", j)
def _fug_coeff_vap(b, j):
return b._fug_coeff_cubic("Vap", j)
def _energy_internal_mol_vap(b):
return b._energy_internal_mol_cubic("Vap")
def _enth_mol_vap(b):
return b._enth_mol_cubic("Vap")
def _enth_mol_vap_ig(b):
return b._enth_mol_ig("Vap")
def _entr_mol_vap(b):
return b._entr_mol_cubic("Vap")
def _entr_mol_vap_ig(b):
return b._entr_mol_ig("Vap")
def bubble_temp_vap(b, j):
def a(k):
return (
b.omegaA
* (
(const.gas_constant * b.params.temperature_crit[k]) ** 2
/ b.params.pressure_crit[k]
)
* (
(
1
+ b.fw[k]
* (
1
- sqrt(b.temperature_bubble / b.params.temperature_crit[k])
)
)
** 2
)
)
am = sum(
sum(
b._mole_frac_tbub[i]
* b._mole_frac_tbub[j]
* sqrt(a(i) * a(j))
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b._mole_frac_tbub[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure / (const.gas_constant * b.temperature_bubble) ** 2
B = bm * b.pressure / (const.gas_constant * b.temperature_bubble)
delta = (
2
* sqrt(a(j))
/ am
* sum(
b._mole_frac_tbub[i] * sqrt(a(i)) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_vap(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
def dew_temp_vap(b, j):
def a(k):
return (
b.omegaA
* (
(const.gas_constant * b.params.temperature_crit[k]) ** 2
/ b.params.pressure_crit[k]
)
* (
(
1
+ b.fw[k]
* (1 - sqrt(b.temperature_dew / b.params.temperature_crit[k]))
)
** 2
)
)
am = sum(
sum(
b.mole_frac_comp[i]
* b.mole_frac_comp[j]
* sqrt(a(i) * a(j))
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b.mole_frac_comp[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure / (const.gas_constant * b.temperature_dew) ** 2
B = bm * b.pressure / (const.gas_constant * b.temperature_dew)
delta = (
2
* sqrt(a(j))
/ am
* sum(
b.mole_frac_comp[i] * sqrt(a(i)) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_vap(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
def bubble_pres_vap(b, j):
am = sum(
sum(
b._mole_frac_pbub[i]
* b._mole_frac_pbub[j]
* sqrt(b.a[i] * b.a[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b._mole_frac_pbub[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure_bubble / (const.gas_constant * b.temperature) ** 2
B = bm * b.pressure_bubble / (const.gas_constant * b.temperature)
delta = (
2
* sqrt(b.a[j])
/ am
* sum(
b._mole_frac_pbub[i] * sqrt(b.a[i]) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_vap(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
def dew_press_vap(b, j):
am = sum(
sum(
b.mole_frac_comp[i]
* b.mole_frac_comp[j]
* sqrt(b.a[i] * b.a[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
bm = sum(b.mole_frac_comp[i] * b.b[i] for i in b.params.component_list)
A = am * b.pressure_dew / (const.gas_constant * b.temperature) ** 2
B = bm * b.pressure_dew / (const.gas_constant * b.temperature)
delta = (
2
* sqrt(b.a[j])
/ am
* sum(
b.mole_frac_comp[i] * sqrt(b.a[i]) * (1 - b.params.kappa[j, i])
for i in b.params.component_list
)
)
Z = b.expr_write.z_vap(eos=b.params.cubic_type, A=A, B=B)
return exp(
(
b.b[j] / bm * (Z - 1) * (B * b.EoS_p)
- safe_log(Z - B, eps=1e-6) * (B * b.EoS_p)
+ A
* (b.b[j] / bm - delta)
* safe_log(
(2 * Z + B * (b.EoS_u + b.EoS_p))
/ (2 * Z + B * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (B * b.EoS_p)
)
# -----------------------------------------------------------------------------
# Common Cubic Functions
# All of these equations drawn from Properties of Gases and Liquids
# Quantities appended with _eq represent calculations @ equilibrium temperature
def common_cubic(blk):
if hasattr(blk, "omegaA"):
return
blk.omegaA = EoS_param[blk.params.cubic_type]["omegaA"]
blk.EoS_Bc = EoS_param[blk.params.cubic_type]["coeff_b"]
blk.EoS_u = EoS_param[blk.params.cubic_type]["u"]
blk.EoS_w = EoS_param[blk.params.cubic_type]["w"]
blk.EoS_p = sqrt(blk.EoS_u**2 - 4 * blk.EoS_w)
# Create expressions for coefficients
def func_fw(b, j):
if b.params.cubic_type == CubicEoS.PR:
return (
0.37464
+ 1.54226 * b.params.omega[j]
- 0.26992 * b.params.omega[j] ** 2
)
elif b.params.cubic_type == CubicEoS.SRK:
return 0.48 + 1.574 * b.params.omega[j] - 0.176 * b.params.omega[j] ** 2
else:
raise BurntToast(
"{} received unrecognised cubic type. This should "
"never happen, so please contact the IDAES developers "
"with this bug.".format(b.name)
)
blk.fw = Param(
blk.params.component_list, initialize=func_fw, doc="EoS S factor"
)
def func_b(b, j):
return (
b.EoS_Bc
* const.gas_constant
* b.params.temperature_crit[j]
/ b.params.pressure_crit[j]
)
blk.b = Param(
blk.params.component_list,
initialize=func_b,
doc="Component b coefficient",
units=pyunits.m**3 / pyunits.mol,
)
def func_a(b, j):
return (
b.omegaA
* (
(const.gas_constant * b.params.temperature_crit[j]) ** 2
/ b.params.pressure_crit[j]
)
* (
(
1
+ b.fw[j]
* (1 - sqrt(b.temperature / b.params.temperature_crit[j]))
)
** 2
)
)
blk.a = Expression(
blk.params.component_list, rule=func_a, doc="Component a coefficient"
)
def func_a_eq(b, j):
return (
b.omegaA
* (
(const.gas_constant * b.params.temperature_crit[j]) ** 2
/ b.params.pressure_crit[j]
)
* (
(1 + b.fw[j] * (1 - sqrt(b._teq / b.params.temperature_crit[j])))
** 2
)
)
blk._a_eq = Expression(
blk.params.component_list,
rule=func_a_eq,
doc="Component a coefficient at Teq",
)
def rule_am(b, p):
return sum(
sum(
b.mole_frac_phase_comp[p, i]
* b.mole_frac_phase_comp[p, j]
* sqrt(b.a[i] * b.a[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
blk.am = Expression(blk.params.phase_list, rule=rule_am)
def rule_am_eq(b, p):
return sum(
sum(
b.mole_frac_phase_comp[p, i]
* b.mole_frac_phase_comp[p, j]
* sqrt(b._a_eq[i] * b._a_eq[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
for i in b.params.component_list
)
blk._am_eq = Expression(blk.params.phase_list, rule=rule_am_eq)
def rule_bm(b, p):
return sum(
b.mole_frac_phase_comp[p, i] * b.b[i] for i in b.params.component_list
)
blk.bm = Expression(blk.params.phase_list, rule=rule_bm)
def rule_A(b, p):
return b.am[p] * b.pressure / (const.gas_constant * b.temperature) ** 2
blk.A = Expression(blk.params.phase_list, rule=rule_A)
def rule_B(b, p):
return b.bm[p] * b.pressure / (const.gas_constant * b.temperature)
blk.B = Expression(blk.params.phase_list, rule=rule_B)
def rule_A_eq(b, p):
return b._am_eq[p] * b.pressure / (const.gas_constant * b._teq) ** 2
blk._A_eq = Expression(blk.params.phase_list, rule=rule_A_eq)
def rule_B_eq(b, p):
return b.bm[p] * b.pressure / (const.gas_constant * b._teq)
blk._B_eq = Expression(blk.params.phase_list, rule=rule_B_eq)
def rule_delta(b, p, i):
# See pg. 145 in Properties of Gases and Liquids
return (
2
* sqrt(blk.a[i])
/ b.am[p]
* sum(
b.mole_frac_phase_comp[p, j]
* sqrt(blk.a[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
)
blk.delta = Expression(
blk.params.phase_list, blk.params.component_list, rule=rule_delta
)
def rule_delta_eq(b, p, i):
# See pg. 145 in Properties of Gases and Liquids
return (
2
* sqrt(blk._a_eq[i])
/ b._am_eq[p]
* sum(
b.mole_frac_phase_comp[p, j]
* sqrt(blk._a_eq[j])
* (1 - b.params.kappa[i, j])
for j in b.params.component_list
)
)
blk._delta_eq = Expression(
blk.params.phase_list, blk.params.component_list, rule=rule_delta_eq
)
def rule_dadT(b, p):
# See pg. 102 in Properties of Gases and Liquids
return -(
(const.gas_constant / 2)
* sqrt(b.omegaA)
* sum(
sum(
b.mole_frac_phase_comp[p, i]
* b.mole_frac_phase_comp[p, j]
* (1 - b.params.kappa[i, j])
* (
b.fw[j]
* sqrt(
b.a[i]
* b.params.temperature_crit[j]
/ b.params.pressure_crit[j]
)
+ b.fw[i]
* sqrt(
b.a[j]
* b.params.temperature_crit[i]
/ b.params.pressure_crit[i]
)
)
for j in b.params.component_list
)
for i in b.params.component_list
)
/ sqrt(b.temperature)
)
blk.dadT = Expression(blk.params.phase_list, rule=rule_dadT)
def rule_compress_fact(b, p):
if p == "Vap":
return blk.expr_write.z_vap(eos=b.params.cubic_type, A=b.A[p], B=b.B[p])
else:
return blk.expr_write.z_liq(eos=b.params.cubic_type, A=b.A[p], B=b.B[p])
blk.compress_fact_phase = Expression(
blk.params.phase_list, rule=rule_compress_fact
)
def rule_compress_fact_eq(b, p):
if p == "Vap":
return blk.expr_write.z_vap(
eos=b.params.cubic_type, A=b._A_eq[p], B=b._B_eq[p]
)
else:
return blk.expr_write.z_liq(
eos=b.params.cubic_type, A=b._A_eq[p], B=b._B_eq[p]
)
blk._compress_fact_eq = Expression(
blk.params.phase_list, rule=rule_compress_fact_eq
)
def _vol_mol_cubic(b, p):
return (
b.pressure * b.vol_mol_phase[p]
== b.compress_fact_phase[p] * const.gas_constant * b.temperature
)
def _dens_mol_cubic(b, p):
return b.pressure == (
b.dens_mol_phase[p]
* b.compress_fact_phase[p]
* const.gas_constant
* b.temperature
)
def _dens_mass_cubic(b, p):
return b.dens_mass_phase[p] == b.dens_mol_phase[p] * b.mw_phase[p]
def _fug_cubic(b, p, j):
return b.mole_frac_phase_comp[p, j] * b.pressure * b.fug_coeff_phase_comp[p, j]
def _fug_coeff_cubic(b, p, j):
# See pg. 145 in Properties of Gases and Liquids
return exp(
(
b.b[j] / b.bm[p] * (b.compress_fact_phase[p] - 1) * (b.B[p] * b.EoS_p)
- safe_log(b.compress_fact_phase[p] - b.B[p], eps=1e-6)
* (b.B[p] * b.EoS_p)
+ b.A[p]
* (b.b[j] / b.bm[p] - b.delta[p, j])
* safe_log(
(2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u + b.EoS_p))
/ (2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
)
/ (b.B[p] * b.EoS_p)
)
def _log_equilibrium_cubic(b, p, j):
# See pg. 145 in Properties of Gases and Liquids
return (
b.b[j] / b.bm[p] * (b._compress_fact_eq[p] - 1) * (b._B_eq[p] * b.EoS_p)
- safe_log(b._compress_fact_eq[p] - b._B_eq[p], eps=1e-6)
* (b._B_eq[p] * b.EoS_p)
+ b._A_eq[p]
* (b.b[j] / b.bm[p] - b._delta_eq[p, j])
* safe_log(
(2 * b._compress_fact_eq[p] + b._B_eq[p] * (b.EoS_u + b.EoS_p))
/ (2 * b._compress_fact_eq[p] + b._B_eq[p] * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
) / (b._B_eq[p] * b.EoS_p) + log(b.mole_frac_phase_comp[p, j])
def _energy_internal_mol_cubic(b, p):
# Derived from equation on pg. 120 in Properties of Gases and Liquids
return (
(b.temperature * b.dadT[p] - b.am[p])
* safe_log(
(2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u + b.EoS_p))
/ (2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
) / (b.bm[p] * b.EoS_p) + b._energy_internal_mol_ig(p)
def _energy_internal_mol_ig(b, p):
dU_form = {}
for j in b.params.component_list:
ele_comp = b.params.get_component(j).config.elemental_composition
if ele_comp is None:
raise ConfigurationError(
"{} calculation of internal energy requires elemental "
"composition of all species. Please set this using the "
"elemental_composition argument in the component "
"declaration ({}).".format(b.name, j)
)
delta_n = 0
for e, s in ele_comp.items():
# Check for any element which is vapor at standard state
if e in ["He", "Ne", "Ar", "Kr", "Xe", "Ra"]:
delta_n += -s
elif e in ["F", "Cl", "H", "N", "O"]:
delta_n += -s / 2 # These are diatomic at standard state
if p == "Vap":
delta_n += 1 # One mole of gaseous compound is formed
dU_form[j] = delta_n * const.gas_constant * b.params.temperature_ref
return sum(
b.mole_frac_phase_comp[p, j]
* (
b._enth_mol_comp_ig(j)
- const.gas_constant * (b.temperature - b.params.temperature_ref)
+ b.params.enth_mol_form_ref[j]
+ dU_form[j]
)
for j in b.params.component_list
)
def _enth_mol_cubic(b, p):
# Derived from equation on pg. 120 in Properties of Gases and Liquids
return (
(b.temperature * b.dadT[p] - b.am[p])
* safe_log(
(2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u + b.EoS_p))
/ (2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
+ const.gas_constant
* b.temperature
* (b.compress_fact_phase[p] - 1)
* b.bm[p]
* b.EoS_p
) / (b.bm[p] * b.EoS_p) + b._enth_mol_ig(p)
def _enth_mol_ig(b, p):
return sum(
b.mole_frac_phase_comp[p, j]
* (b._enth_mol_comp_ig(j) + b.params.enth_mol_form_ref[j])
for j in b.params.component_list
)
def _entr_mol_cubic(b, p):
# See pg. 102 in Properties of Gases and Liquids
return (
const.gas_constant
* safe_log(
(b.compress_fact_phase[p] - b.B[p]) / b.compress_fact_phase[p], eps=1e-6
)
* b.bm[p]
* b.EoS_p
+ const.gas_constant
* safe_log(
b.compress_fact_phase[p] * b.params.pressure_ref / b.pressure, eps=1e-6
)
* b.bm[p]
* b.EoS_p
+ b.dadT[p]
* safe_log(
(2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u + b.EoS_p))
/ (2 * b.compress_fact_phase[p] + b.B[p] * (b.EoS_u - b.EoS_p)),
eps=1e-6,
)
) / (b.bm[p] * b.EoS_p) + b._entr_mol_ig(p)
def _entr_mol_ig(b, p):
return sum(
b.mole_frac_phase_comp[p, j]
* (b._entr_mol_comp_ig(j) + b.params.entr_mol_form_ref[j])
for j in b.params.component_list
)
# -----------------------------------------------------------------------------
# Pure component properties
def _enth_mol_comp_ig(b, j):
return (
(b.params.cp_mol_ig_comp_coeff_4[j] / 4)
* (b.temperature**4 - b.params.temperature_ref**4)
+ (b.params.cp_mol_ig_comp_coeff_3[j] / 3)
* (b.temperature**3 - b.params.temperature_ref**3)
+ (b.params.cp_mol_ig_comp_coeff_2[j] / 2)
* (b.temperature**2 - b.params.temperature_ref**2)
+ b.params.cp_mol_ig_comp_coeff_1[j]
* (b.temperature - b.params.temperature_ref)
)
def _entr_mol_comp_ig(b, j):
return (
(b.params.cp_mol_ig_comp_coeff_4[j] / 3)
* (b.temperature**3 - b.params.temperature_ref**3)
+ (b.params.cp_mol_ig_comp_coeff_3[j] / 2)
* (b.temperature**2 - b.params.temperature_ref**2)
+ b.params.cp_mol_ig_comp_coeff_2[j]
* (b.temperature - b.params.temperature_ref)
+ b.params.cp_mol_ig_comp_coeff_1[j]
* log(b.temperature / b.params.temperature_ref)
)
def calculate_scaling_factors(self):
# Get default scale factors and do calculations from base classes
super().calculate_scaling_factors()
phases = self.params.config.valid_phase
is_two_phase = len(phases) == 2 # possibly {Liq}, {Vap}, or {Liq, Vap}
sf_flow = iscale.get_scaling_factor(self.flow_mol, default=1, warning=True)
sf_T = iscale.get_scaling_factor(self.temperature, default=1, warning=True)
if self.is_property_constructed("_teq"):
iscale.set_scaling_factor(self._teq, sf_T)
if self.is_property_constructed("_teq_constraint"):
iscale.constraint_scaling_transform(
self._teq_constraint, sf_T, overwrite=False
)
if self.is_property_constructed("_t1"):
iscale.set_scaling_factor(self._t1, sf_T)
if self.is_property_constructed("_t1_constraint"):
iscale.constraint_scaling_transform(
self._t1_constraint, sf_T, overwrite=False
)
if self.is_property_constructed("_mole_frac_pdew"):
iscale.set_scaling_factor(self._mole_frac_pdew, 1e3)
iscale.constraint_scaling_transform(
self._sum_mole_frac_pdew, 1e3, overwrite=False
)
if self.is_property_constructed("total_flow_balance"):
s = iscale.get_scaling_factor(self.flow_mol, default=1, warning=True)
iscale.constraint_scaling_transform(
self.total_flow_balance, s, overwrite=False
)
if self.is_property_constructed("component_flow_balances"):
for i, c in self.component_flow_balances.items():
if is_two_phase:
s = iscale.get_scaling_factor(
self.mole_frac_comp[i], default=1, warning=True
)
s *= sf_flow
iscale.constraint_scaling_transform(c, s, overwrite=False)
else:
s = iscale.get_scaling_factor(
self.mole_frac_comp[i], default=1, warning=True
)
iscale.constraint_scaling_transform(c, s, overwrite=False)
if self.is_property_constructed("dens_mol_phase"):
for c in self.eq_dens_mol_phase.values():
sf = iscale.get_scaling_factor(
self.dens_mol_phase, default=1, warning=True
)
iscale.constraint_scaling_transform(c, sf, overwrite=False)
if self.is_property_constructed("dens_mass_phase"):
for p, c in self.eq_dens_mass_phase.items():
sf = iscale.get_scaling_factor(
self.dens_mass_phase[p], default=1, warning=True
)
iscale.constraint_scaling_transform(c, sf, overwrite=False)
if self.is_property_constructed("enth_mol_phase"):
for p, c in self.eq_enth_mol_phase.items():
sf = iscale.get_scaling_factor(
self.enth_mol_phase[p], default=1, warning=True
)
iscale.constraint_scaling_transform(c, sf, overwrite=False)
if self.is_property_constructed("enth_mol"):
sf = iscale.get_scaling_factor(self.enth_mol, default=1, warning=True)
sf *= sf_flow
iscale.constraint_scaling_transform(self.eq_enth_mol, sf, overwrite=False)
if self.is_property_constructed("entr_mol_phase"):
for p, c in self.eq_entr_mol_phase.items():
sf = iscale.get_scaling_factor(
self.entr_mol_phase[p], default=1, warning=True
)
iscale.constraint_scaling_transform(c, sf, overwrite=False)
if self.is_property_constructed("entr_mol"):
sf = iscale.get_scaling_factor(self.entr_mol, default=1, warning=True)
sf *= sf_flow
iscale.constraint_scaling_transform(self.eq_entr_mol, sf, overwrite=False)
if self.is_property_constructed("gibbs_mol_phase"):
for p, c in self.eq_gibbs_mol_phase.items():
sf = iscale.get_scaling_factor(
self.gibbs_mol_phase[p], default=1, warning=True
)
iscale.constraint_scaling_transform(c, sf, overwrite=False)
