diff --git a/idaes/core/util/initialization.py b/idaes/core/util/initialization.py index 3436d2c5ea..4ddca90536 100644 --- a/idaes/core/util/initialization.py +++ b/idaes/core/util/initialization.py @@ -21,6 +21,7 @@ Constraint, value, ) +from pyomo.common.collections import ComponentMap from pyomo.network import Arc from pyomo.dae import ContinuousSet from pyomo.core.expr.visitor import identify_variables @@ -533,3 +534,35 @@ def initialize_by_time_element(fs, time, **kwargs): # Logger message that initialization is finished init_log.info("Initialization completed. Model has been reactivated") + +def _fix_vars(var_list): + flags = ComponentMap() + for var in var_list: + if var.is_indexed(): + for subvar in var.values(): + flags[subvar] = subvar.fixed + subvar.fix() + else: + flags[var] = var.fixed + var.fix() + return flags + + +def _unfix_vars(var_list): + flags = ComponentMap() + if var.is_indexed(): + for subvar in var.values(): + flags[subvar] = subvar.fixed + subvar.unfix() + else: + flags[var] = var.fixed + var.unfix() + return flags + + +def _restore_fixedness(flags): + for var, fix in flags.items(): + if fix: + var.fix() + else: + var.unfix() \ No newline at end of file diff --git a/idaes/models/properties/modular_properties/base/generic_property.py b/idaes/models/properties/modular_properties/base/generic_property.py index f90609b0fb..81366f11e4 100644 --- a/idaes/models/properties/modular_properties/base/generic_property.py +++ b/idaes/models/properties/modular_properties/base/generic_property.py @@ -4689,7 +4689,7 @@ def _log_mole_frac_phase_comp(self): try: self.log_mole_frac_phase_comp = Var( self.phase_component_set, - initialize=0, + initialize=-1, bounds=(-100, log(1.001)), units=pyunits.dimensionless, doc="Log of mole fractions of component by phase", @@ -4714,8 +4714,8 @@ def _log_mole_frac_phase_comp_apparent(self): try: self.log_mole_frac_phase_comp_apparent = Var( self.params.apparent_phase_component_set, - initialize=0, - bounds=(-50, 0), + initialize=-1, + bounds=(-100, log(1.001)), units=pyunits.dimensionless, doc="Log of mole fractions of component by phase", ) @@ -4739,8 +4739,8 @@ def _log_mole_frac_phase_comp_true(self): try: self.log_mole_frac_phase_comp_true = Var( self.params.true_phase_component_set, - initialize=0, - bounds=(-50, 0), + initialize=-1, + bounds=(-100, log(1.001)), units=pyunits.dimensionless, doc="Log of mole fractions of component by phase", ) diff --git a/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py b/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py index 391a8e92f9..a1acadbec4 100644 --- a/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py +++ b/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py @@ -633,8 +633,10 @@ def test_equilibrium_form(self, model): ) assert value(model.rblock[1].equilibrium_constraint["e2"].body) == value( model.rblock[1].log_k_eq["e2"] - - model.sblock[1].log_mole_frac_phase_comp["p2", "c1"] * -5 - + model.sblock[1].log_mole_frac_phase_comp["p2", "c2"] * 6 + - ( + model.sblock[1].log_mole_frac_phase_comp["p2", "c1"] * -5 + + model.sblock[1].log_mole_frac_phase_comp["p2", "c2"] * 6 + ) ) @pytest.mark.unit diff --git a/idaes/models/properties/modular_properties/state_definitions/electrolyte_states.py b/idaes/models/properties/modular_properties/state_definitions/electrolyte_states.py index cdb048a433..311219856c 100644 --- a/idaes/models/properties/modular_properties/state_definitions/electrolyte_states.py +++ b/idaes/models/properties/modular_properties/state_definitions/electrolyte_states.py @@ -118,7 +118,7 @@ def appr_to_true_species(b, p, j): # Next, check for inherent reactions if b.has_inherent_reactions: for r in b.params.inherent_reaction_idx: - # Get stoichiometric coeffiicient for inherent reactions + # Get stoichiometric coefficient for inherent reactions gamma = b.params.inherent_reaction_stoichiometry[r, p, j] if gamma != 0: @@ -156,16 +156,37 @@ def true_species_mole_fractions(b, p, j): def _apparent_species_scaling(b): for p, j in b.params.true_phase_component_set: - sf = iscale.get_scaling_factor( - b.flow_mol_phase_comp_true[p, j], default=1, warning=True - ) + pobj = b.params.get_phase(p) + cobj = b.params.get_component(j) - iscale.constraint_scaling_transform( - b.appr_to_true_species[p, j], sf, overwrite=False - ) - iscale.constraint_scaling_transform( - b.true_mole_frac_constraint[p, j], sf, overwrite=False - ) + # For apparent species in an aqueous phase, want to scale material balance + # equations by apparent magnitude. If they're not consumed to exhaustion, + # they'll be well-scaled in this equation. If they are consumed to "exhaustion" + # (several orders of magnitude below apparent concentration), then their concentration + # will decouple from this equation and will be given by equilibrium equation instead + if pobj.is_aqueous_phase and not isinstance(cobj, IonData): + sf_apparent = iscale.get_scaling_factor( + b.flow_mol_phase_comp_apparent[p, j], default=1, warning=True + ) + iscale.constraint_scaling_transform( + b.appr_to_true_species[p, j], sf_apparent, overwrite=False + ) + sf_true = iscale.get_scaling_factor( + b.flow_mol_phase_comp_true[p, j], default=1, warning=True + ) + iscale.constraint_scaling_transform( + b.true_mole_frac_constraint[p, j], sf_true, overwrite=False + ) + else: + sf = iscale.get_scaling_factor( + b.flow_mol_phase_comp_true[p, j], default=1, warning=True + ) + iscale.constraint_scaling_transform( + b.appr_to_true_species[p, j], sf, overwrite=False + ) + iscale.constraint_scaling_transform( + b.true_mole_frac_constraint[p, j], sf, overwrite=False + ) def _true_species_state(b): diff --git a/idaes/models/unit_models/heat_exchanger_ntu.py b/idaes/models/unit_models/heat_exchanger_ntu.py index 956af7b52c..a05ec6ef22 100644 --- a/idaes/models/unit_models/heat_exchanger_ntu.py +++ b/idaes/models/unit_models/heat_exchanger_ntu.py @@ -48,6 +48,7 @@ from idaes.core.util.math import smooth_min, smooth_max from idaes.core.solvers import get_solver from idaes.core.util.exceptions import InitializationError +import idaes.core.util.scaling as iscale import idaes.logger as idaeslog from idaes.core.initialization import SingleControlVolumeUnitInitializer @@ -478,7 +479,13 @@ def rule_entu(blk, t): self.heat_duty_constraint = Constraint(self.flowsheet().time, rule=rule_entu) - # TODO : Add scaling methods + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + + for t in self.flowsheet().time: + sf_heat = iscale.get_scaling_factor(self.hot_side.heat[t], default=1, warning=True) + iscale.constraint_scaling_transform(self.energy_balance_constraint[t], sf_heat, overwrite=False) + iscale.constraint_scaling_transform(self.heat_duty_constraint[t], sf_heat, overwrite=False) def initialize_build( self, diff --git a/idaes/models/unit_models/stream_scaler.py b/idaes/models/unit_models/stream_scaler.py new file mode 100644 index 0000000000..387a6f6641 --- /dev/null +++ b/idaes/models/unit_models/stream_scaler.py @@ -0,0 +1,299 @@ +################################################################################# +# 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-2023 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. +################################################################################# +""" +Unit model to adjust size of streams to represent, for example, a stream being split across several identical units, +which are then all modeled as a single IDAES unit +""" +from enum import Enum +from functools import partial + +from pyomo.environ import ( + Block, + check_optimal_termination, + Param, + PositiveReals, + Reals, + RangeSet, + units as pyunits, + Var, +) +from pyomo.network import Port +from pyomo.common.config import ConfigBlock, ConfigValue, In, ListOf, Bool + +from idaes.core import ( + declare_process_block_class, + UnitModelBlockData, + useDefault, + MaterialBalanceType, + MaterialFlowBasis, +) +from idaes.core.util.config import ( + is_physical_parameter_block, + is_state_block, +) +from idaes.core.util.exceptions import ( + BurntToast, + ConfigurationError, + PropertyNotSupportedError, + InitializationError, +) +from idaes.core.base.var_like_expression import VarLikeExpression +from idaes.core.util.math import smooth_min +from idaes.core.util.tables import create_stream_table_dataframe +import idaes.core.util.scaling as iscale +from idaes.core.solvers import get_solver +import idaes.logger as idaeslog + +__author__ = "Douglas Allan" + + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("StreamScaler") +class StreamScalerData(UnitModelBlockData): + """ + This is a general purpose model for a Mixer block with the IDAES modeling + framework. This block can be used either as a stand-alone Mixer unit + operation, or as a sub-model within another unit operation. + + This model creates a number of StateBlocks to represent the incoming + streams, then writes a set of phase-component material balances, an + overall enthalpy balance and a momentum balance (2 options) linked to a + mixed-state StateBlock. The mixed-state StateBlock can either be specified + by the user (allowing use as a sub-model), or created by the Mixer. + + When being used as a sub-model, Mixer should only be used when a set + of new StateBlocks are required for the streams to be mixed. It should not + be used to mix streams from mutiple ControlVolumes in a single unit model - + in these cases the unit model developer should write their own mixing + equations. + """ + + CONFIG = ConfigBlock() + CONFIG.declare( + "dynamic", + ConfigValue( + domain=In([False]), + default=False, + description="Dynamic model flag - must be False", + doc="""Indicates whether this model will be dynamic or not, +**default** = False. Scalar blocks are always steady-state.""", + ), + ) + CONFIG.declare( + "has_holdup", + ConfigValue( + default=False, + domain=In([False]), + description="Holdup construction flag - must be False", + doc="""Scalar blocks do not contain holdup, thus this must be +False.""", + ), + ) + CONFIG.declare( + "property_package", + ConfigValue( + default=useDefault, + domain=is_physical_parameter_block, + description="Property package to use for mixer", + doc="""Property parameter object used to define property +calculations, **default** - useDefault. +**Valid values:** { +**useDefault** - use default package from parent model or flowsheet, +**PropertyParameterObject** - a PropertyParameterBlock object.}""", + ), + ) + CONFIG.declare( + "property_package_args", + ConfigBlock( + implicit=True, + description="Arguments to use for constructing property packages", + doc="""A ConfigBlock with arguments to be passed to a property +block(s) and used when constructing these, +**default** - None. +**Valid values:** { +see property package for documentation.}""", + ), + ) +# CONFIG.declare( +# "mixed_state_block", +# ConfigValue( +# default=None, +# domain=is_state_block, +# description="Existing StateBlock to use as mixed stream", +# doc="""An existing state block to use as the outlet stream from the +# Mixer block, +# **default** - None. +# **Valid values:** { +# **None** - create a new StateBlock for the mixed stream, +# **StateBlock** - a StateBock to use as the destination for the mixed stream.} +# """, +# ), +# ) +# CONFIG.declare( +# "construct_ports", +# ConfigValue( +# default=True, +# domain=Bool, +# description="Construct inlet and outlet Port objects", +# doc="""Argument indicating whether model should construct Port +# objects linked to all inlet states and the mixed state, +# **default** - True. +# **Valid values:** { +# **True** - construct Ports for all states, +# **False** - do not construct Ports.""", +# ), +# ) + + def build(self): + """ + General build method for StreamScalarData. This method calls a number + of sub-methods which automate the construction of expected attributes + of unit models. + + Inheriting models should call `super().build`. + + Args: + None + + Returns: + None + """ + # Call super.build() + super(StreamScalerData, self).build() + + tmp_dict = dict(**self.config.property_package_args) + tmp_dict["has_phase_equilibrium"] = False + tmp_dict["defined_state"] = True + + # Call setup methods from ControlVolumeBlockData + self._get_property_package() + self._get_indexing_sets() + + self.inlet_block = self.config.property_package.build_state_block( + self.flowsheet().time, doc="Material properties at inlet", **tmp_dict + ) + self.outlet_block = Block() + self.multiplier = Var( + initialize=1, + domain=PositiveReals, + units=pyunits.dimensionless, + doc="Factor by which to scale dimensionless streams" + ) + self.add_inlet_port(name="inlet", block=self.inlet_block) + self.outlet = Port(doc="Outlet port") + + def rule_scale_var(b, *args, var=None): + return self.multiplier * var[args] + + def rule_no_scale_var(b, *args, var=None): + return var[args] + + for var_name in self.inlet.vars.keys(): + var = getattr(self.inlet, var_name) + if "flow" in var_name: + rule=partial(rule_scale_var, var=var) + else: + rule=partial(rule_no_scale_var, var=var) + self.outlet_block.add_component( + var_name, + VarLikeExpression( + var.index_set(), + rule=rule + ) + ) + expr = getattr(self.outlet_block, var_name) + self.outlet.add(expr, var_name) + + def initialize_build( + blk, outlvl=idaeslog.NOTSET, optarg=None, solver=None, hold_state=False + ): + """ + Initialization routine for mixer. + + Keyword Arguments: + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None, use + default solver options) + 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. **Valid values:** + **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="unit") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") + + # Create solver + + # Initialize inlet state blocks + flags = blk.inlet_block.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + hold_state=True, + ) + + if hold_state is True: + return flags + else: + blk.release_state(flags, outlvl=outlvl) + + 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 logging + + Returns: + None + """ + blk.inlet_block.release_state(flags, outlvl=outlvl) + + def _get_stream_table_contents(self, time_point=0): + io_dict = { + "Inlet": self.inlet, + "Outlet": self.outlet, + } + return create_stream_table_dataframe(io_dict, time_point=time_point) + + def calculate_scaling_factors(self): + # Scaling factors for the property block are calculated automatically + super().calculate_scaling_factors() + + # Need to pass on scaling factors from the property block to the outlet + # VarLikeExpressions so arcs get scaled right + scale = 1/self.multiplier.value + for var_name in self.inlet.vars.keys(): + var = getattr(self.inlet, var_name) + outlet_expr = getattr(self.outlet, var_name) + for key, subvar in var.items(): + sf = iscale.get_scaling_factor(subvar, default=1, warning=True) + iscale.set_scaling_factor(outlet_expr[key],scale*sf) + \ No newline at end of file diff --git a/idaes/models/unit_models/tests/test_stream_scaler.py b/idaes/models/unit_models/tests/test_stream_scaler.py new file mode 100644 index 0000000000..a25388b104 --- /dev/null +++ b/idaes/models/unit_models/tests/test_stream_scaler.py @@ -0,0 +1,1035 @@ +################################################################################# +# 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-2023 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. +################################################################################# +""" +Tests for StreamScaler. + +Author: Douglas Allan, Andrew Lee +""" +import pytest +import pandas + +from pyomo.environ import ( + Block, + check_optimal_termination, + ConcreteModel, + Constraint, + Param, + RangeSet, + Set, + Var, + value, + units as pyunits, +) +from pyomo.util.check_units import assert_units_consistent + +from pyomo.network import Port +from pyomo.common.config import ConfigBlock + +from idaes.core import ( + FlowsheetBlock, + declare_process_block_class, + StateBlockData, + StateBlock, + PhysicalParameterBlock, + MaterialBalanceType, + EnergyBalanceType, + Phase, + Component, +) +from idaes.models.properties.activity_coeff_models.BTX_activity_coeff_VLE import ( + BTXParameterBlock, +) +from idaes.models.properties import iapws95 +from idaes.models.properties.examples.saponification_thermo import ( + SaponificationParameterBlock, +) + +from idaes.models.unit_models.stream_scaler import ( + StreamScaler, + StreamScalerData, + MixingType, + MomentumMixingType, +) +from idaes.core.util.exceptions import ( + BurntToast, + ConfigurationError, + InitializationError, + PropertyNotSupportedError, +) +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_variables, + number_total_constraints, + number_unused_variables, +) +from idaes.core.util.testing import ( + PhysicalParameterTestBlock, + TestStateBlock, + initialization_tester, +) +from idaes.models.properties.modular_properties.base.generic_property import ( + GenericParameterBlock, +) +from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available +from idaes.core.base.var_like_expression import VarLikeExpression + +# TODO: Should have a test for this that does not requrie models_extra +from idaes.models_extra.power_generation.properties.natural_gas_PR import get_prop +import idaes.core.util.scaling as iscale +from idaes.core.solvers import get_solver + + +# ----------------------------------------------------------------------------- +# Get default solver for testing +solver = get_solver() + + +# ----------------------------------------------------------------------------- +# Unit Tests for StreamScaler +class TestStreamScaler(object): + @declare_process_block_class("StreamScalerFrame") + class StreamScalerFrameData(StreamScalerData): + def build(self): + super(StreamScalerData, self).build() + + @pytest.fixture(scope="function") + def scaler_frame(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.scaler = StreamScalerFrame(property_package=m.fs.pp) + + return m + + @pytest.mark.unit + def test_scaler_config(self, scaler_frame): + assert len(scaler_frame.fs.scaler.config) == 4 + assert scaler_frame.fs.scaler.config.dynamic is False + assert scaler_frame.fs.scaler.config.has_holdup is False + assert scaler_frame.fs.scaler.config.property_package == scaler_frame.fs.pp + assert isinstance(scaler_frame.fs.scaler.config.property_package_args, ConfigBlock) + assert len(scaler_frame.fs.scaler.config.property_package_args) == 0 + + @pytest.mark.unit + def test_inherited_methods(self, scaler_frame): + scaler_frame.fs.scaler._get_property_package() + scaler_frame.fs.scaler._get_indexing_sets() + + assert hasattr(scaler_frame.fs.scaler.config.property_package, "phase_list") + + @pytest.mark.unit + def test_add_port_objects(self, scaler_frame): + scaler_frame.fs.scaler._get_property_package() + scaler_frame.fs.scaler._get_indexing_sets() + + inlet_list = scaler_frame.fs.scaler.create_inlet_list() + inlet_blocks = scaler_frame.fs.scaler.add_inlet_state_blocks(inlet_list) + mixed_block = scaler_frame.fs.scaler.add_mixed_state_block() + + scaler_frame.fs.scaler.add_port_objects(inlet_list, inlet_blocks, mixed_block) + + assert isinstance(scaler_frame.fs.scaler.inlet_1, Port) + assert isinstance(scaler_frame.fs.scaler.inlet_2, Port) + assert isinstance(scaler_frame.fs.scaler.outlet, Port) + + # ------------------------------------------------------------------------- + # Test build method + @pytest.mark.build + @pytest.mark.unit + def test_build_default(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.scaler = StreamScaler(property_package=m.fs.pp) + + assert isinstance(m.fs.scaler.scaling_factor, Var) + assert len(m.fs.scaler.scaling_factor) == 1 + + assert isinstance(m.fs.scaler.outlet_block, Block) + assert len(m.fs.scaler.outlet_block) == 1 + + assert isinstance(m.fs.scaler.inlet_1, Port) + assert isinstance(m.fs.scaler.inlet_2, Port) + assert isinstance(m.fs.scaler.outlet, Port) + + @pytest.mark.build + @pytest.mark.unit + def test_build_phase_equilibrium(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.scaler = Mixer(property_package=m.fs.pp, has_phase_equilibrium=True) + + assert isinstance(m.fs.scaler.material_mixing_equations, Constraint) + assert len(m.fs.scaler.material_mixing_equations) == 4 + assert isinstance(m.fs.scaler.phase_equilibrium_generation, Var) + + assert isinstance(m.fs.scaler.enthalpy_mixing_equations, Constraint) + assert len(m.fs.scaler.enthalpy_mixing_equations) == 1 + + assert m.fs.scaler.inlet_idx.ctype is RangeSet + assert isinstance(m.fs.scaler.minimum_pressure, Var) + assert len(m.fs.scaler.minimum_pressure) == 2 + assert isinstance(m.fs.scaler.eps_pressure, Param) + assert isinstance(m.fs.scaler.minimum_pressure_constraint, Constraint) + assert len(m.fs.scaler.minimum_pressure) == 2 + assert isinstance(m.fs.scaler.mixture_pressure, Constraint) + + assert isinstance(m.fs.scaler.inlet_1, Port) + assert isinstance(m.fs.scaler.inlet_2, Port) + assert isinstance(m.fs.scaler.outlet, Port) + + @pytest.mark.build + @pytest.mark.unit + def test_build_phase_pressure_equality(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.scaler = Mixer( + property_package=m.fs.pp, momentum_mixing_type=MomentumMixingType.equality + ) + + assert isinstance(m.fs.scaler.material_mixing_equations, Constraint) + assert len(m.fs.scaler.material_mixing_equations) == 4 + + assert isinstance(m.fs.scaler.enthalpy_mixing_equations, Constraint) + assert len(m.fs.scaler.enthalpy_mixing_equations) == 1 + + assert isinstance(m.fs.scaler.pressure_equality_constraints, Constraint) + assert len(m.fs.scaler.pressure_equality_constraints) == 2 + + assert isinstance(m.fs.scaler.inlet_1, Port) + assert isinstance(m.fs.scaler.inlet_2, Port) + assert isinstance(m.fs.scaler.outlet, Port) + + # ------------------------------------------------------------------------- + # Test models checks, initialize and release state methods + @pytest.mark.unit + def test_model_checks(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.scaler = Mixer( + property_package=m.fs.pp, momentum_mixing_type=MomentumMixingType.equality + ) + + m.fs.scaler.model_check() + + assert m.fs.scaler.inlet_1_state[0].check is True + assert m.fs.scaler.inlet_2_state[0].check is True + assert m.fs.scaler.mixed_state[0].check is True + + @pytest.mark.ui + @pytest.mark.unit + def test_get_stream_table_contents(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.sb = TestStateBlock(m.fs.time, parameters=m.fs.pp) + + m.fs.scaler = Mixer(property_package=m.fs.pp) + + stable = m.fs.scaler._get_stream_table_contents() + + expected = pandas.DataFrame.from_dict( + { + "Units": { + "component_flow_phase ('p1', 'c1')": getattr( + pyunits.pint_registry, "mole/second" + ), + "component_flow_phase ('p1', 'c2')": getattr( + pyunits.pint_registry, "mole/second" + ), + "component_flow_phase ('p2', 'c1')": getattr( + pyunits.pint_registry, "mole/second" + ), + "component_flow_phase ('p2', 'c2')": getattr( + pyunits.pint_registry, "mole/second" + ), + "temperature": getattr(pyunits.pint_registry, "K"), + "pressure": getattr(pyunits.pint_registry, "Pa"), + }, + "inlet_1": { + "component_flow_phase ('p1', 'c1')": 2.00, + "component_flow_phase ('p1', 'c2')": 2.00, + "component_flow_phase ('p2', 'c1')": 2.00, + "component_flow_phase ('p2', 'c2')": 2.00, + "temperature": 300, + "pressure": 1e5, + }, + "inlet_2": { + "component_flow_phase ('p1', 'c1')": 2.00, + "component_flow_phase ('p1', 'c2')": 2.00, + "component_flow_phase ('p2', 'c1')": 2.00, + "component_flow_phase ('p2', 'c2')": 2.00, + "temperature": 300, + "pressure": 1e5, + }, + "Outlet": { + "component_flow_phase ('p1', 'c1')": 2.00, + "component_flow_phase ('p1', 'c2')": 2.00, + "component_flow_phase ('p2', 'c1')": 2.00, + "component_flow_phase ('p2', 'c2')": 2.00, + "temperature": 300, + "pressure": 1e5, + }, + } + ) + + pandas.testing.assert_frame_equal(stable, expected, rtol=1e-4, atol=1e-4) + + @pytest.mark.initialization + @pytest.mark.component + def test_initialize(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + m.fs.sb = TestStateBlock(m.fs.time, parameters=m.fs.pp) + + m.fs.scaler = Mixer(property_package=m.fs.pp, mixed_state_block=m.fs.sb) + + # Change one inlet pressure to check initialization calculations + m.fs.scaler.inlet_1_state[0].pressure = 8e4 + + f = m.fs.scaler.initialize(hold_state=True) + + assert m.fs.scaler.inlet_1_state[0].init_test is True + assert m.fs.scaler.inlet_2_state[0].init_test is True + assert m.fs.sb[0].init_test is True + assert m.fs.scaler.inlet_1_state[0].hold_state is True + assert m.fs.scaler.inlet_2_state[0].hold_state is True + assert m.fs.sb[0].hold_state is False + + assert m.fs.sb[0].flow_mol_phase_comp["p1", "c1"].value == 4 + assert m.fs.sb[0].flow_mol_phase_comp["p1", "c2"].value == 4 + assert m.fs.sb[0].flow_mol_phase_comp["p2", "c1"].value == 4 + assert m.fs.sb[0].flow_mol_phase_comp["p2", "c2"].value == 4 + + assert m.fs.sb[0].temperature.value == 300 + + assert m.fs.sb[0].pressure.value == 8e4 + + m.fs.scaler.release_state(flags=f) + + assert m.fs.scaler.inlet_1_state[0].hold_state is False + assert m.fs.scaler.inlet_2_state[0].hold_state is False + assert m.fs.sb[0].hold_state is False + + +# ----------------------------------------------------------------------------- +class TestBTX(object): + @pytest.fixture(scope="class") + def btx(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = BTXParameterBlock(valid_phase="Liq") + + m.fs.unit = Mixer(property_package=m.fs.properties) + + m.fs.unit.inlet_1.flow_mol[0].fix(5) # mol/s + m.fs.unit.inlet_1.temperature[0].fix(365) # K + m.fs.unit.inlet_1.pressure[0].fix(2e5) # Pa + m.fs.unit.inlet_1.mole_frac_comp[0, "benzene"].fix(0.5) + m.fs.unit.inlet_1.mole_frac_comp[0, "toluene"].fix(0.5) + + m.fs.unit.inlet_2.flow_mol[0].fix(1) # mol/s + m.fs.unit.inlet_2.temperature[0].fix(300) # K + m.fs.unit.inlet_2.pressure[0].fix(101325) # Pa + m.fs.unit.inlet_2.mole_frac_comp[0, "benzene"].fix(0.5) + m.fs.unit.inlet_2.mole_frac_comp[0, "toluene"].fix(0.5) + + return m + + @pytest.mark.build + @pytest.mark.unit + def test_build(self, btx): + assert hasattr(btx.fs.unit, "inlet_1") + assert len(btx.fs.unit.inlet_1.vars) == 4 + assert hasattr(btx.fs.unit.inlet_1, "flow_mol") + assert hasattr(btx.fs.unit.inlet_1, "mole_frac_comp") + assert hasattr(btx.fs.unit.inlet_1, "temperature") + assert hasattr(btx.fs.unit.inlet_1, "pressure") + + assert hasattr(btx.fs.unit, "inlet_2") + assert len(btx.fs.unit.inlet_2.vars) == 4 + assert hasattr(btx.fs.unit.inlet_2, "flow_mol") + assert hasattr(btx.fs.unit.inlet_2, "mole_frac_comp") + assert hasattr(btx.fs.unit.inlet_2, "temperature") + assert hasattr(btx.fs.unit.inlet_2, "pressure") + + assert hasattr(btx.fs.unit, "outlet") + assert len(btx.fs.unit.outlet.vars) == 4 + assert hasattr(btx.fs.unit.outlet, "flow_mol") + assert hasattr(btx.fs.unit.outlet, "mole_frac_comp") + assert hasattr(btx.fs.unit.outlet, "temperature") + assert hasattr(btx.fs.unit.outlet, "pressure") + + assert number_variables(btx) == 35 + assert number_total_constraints(btx) == 25 + assert number_unused_variables(btx) == 0 + + @pytest.mark.component + def test_units(self, btx): + assert_units_consistent(btx) + + @pytest.mark.unit + def test_dof(self, btx): + assert degrees_of_freedom(btx) == 0 + + @pytest.mark.ui + @pytest.mark.unit + def test_get_stream_table_contents(self, btx): + stable = btx.fs.unit._get_stream_table_contents() + + expected = pandas.DataFrame.from_dict( + { + "Units": { + "flow_mol": getattr(pyunits.pint_registry, "mole/second"), + "mole_frac_comp benzene": getattr( + pyunits.pint_registry, "dimensionless" + ), + "mole_frac_comp toluene": getattr( + pyunits.pint_registry, "dimensionless" + ), + "temperature": getattr(pyunits.pint_registry, "kelvin"), + "pressure": getattr(pyunits.pint_registry, "Pa"), + }, + "inlet_1": { + "flow_mol": 5.0, + "mole_frac_comp benzene": 0.5, + "mole_frac_comp toluene": 0.5, + "temperature": 365, + "pressure": 2e5, + }, + "inlet_2": { + "flow_mol": 1.0, + "mole_frac_comp benzene": 0.5, + "mole_frac_comp toluene": 0.5, + "temperature": 300, + "pressure": 101325.0, + }, + "Outlet": { + "flow_mol": 1.0, + "mole_frac_comp benzene": 0.5, + "mole_frac_comp toluene": 0.5, + "temperature": 298.15, + "pressure": 101325.0, + }, + } + ) + + pandas.testing.assert_frame_equal(stable, expected, rtol=1e-4, atol=1e-4) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_initialize(self, btx): + initialization_tester(btx) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solve(self, btx): + results = solver.solve(btx) + + # Check for optimal solution + assert check_optimal_termination(results) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solution(self, btx): + assert pytest.approx(6, abs=1e-3) == value(btx.fs.unit.outlet.flow_mol[0]) + assert pytest.approx(354.7, abs=1e-1) == value( + btx.fs.unit.outlet.temperature[0] + ) + assert pytest.approx(101325, abs=1e3) == value(btx.fs.unit.outlet.pressure[0]) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_conservation(self, btx): + assert ( + abs( + value( + btx.fs.unit.inlet_1.flow_mol[0] + + btx.fs.unit.inlet_2.flow_mol[0] + - btx.fs.unit.outlet.flow_mol[0] + ) + ) + <= 1e-6 + ) + + assert 1e-6 >= abs( + value( + btx.fs.unit.inlet_1.flow_mol[0] + * btx.fs.unit.inlet_1_state[0].enth_mol_phase["Liq"] + + btx.fs.unit.inlet_2.flow_mol[0] + * btx.fs.unit.inlet_2_state[0].enth_mol_phase["Liq"] + - btx.fs.unit.outlet.flow_mol[0] + * btx.fs.unit.mixed_state[0].enth_mol_phase["Liq"] + ) + ) + + +# ----------------------------------------------------------------------------- +# Tests for Mixer in cases where proeprties do not support pressure +@declare_process_block_class("NoPressureTestBlock") +class _NoPressureParameterBlock(PhysicalParameterBlock): + def build(self): + super(_NoPressureParameterBlock, self).build() + + self.p1 = Phase() + self.p2 = Phase() + self.c1 = Component() + self.c2 = Component() + + self.phase_equilibrium_idx = Set(initialize=["e1", "e2"]) + + self.phase_equilibrium_list = { + "e1": ["c1", ("p1", "p2")], + "e2": ["c2", ("p1", "p2")], + } + + self._state_block_class = NoPressureStateBlock + + @classmethod + def define_metadata(cls, obj): + obj.add_properties({}) + obj.add_default_units( + { + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.g, + "amount": pyunits.mol, + "temperature": pyunits.K, + } + ) + + +@declare_process_block_class("NoPressureStateBlock", block_class=StateBlock) +class NoPressureStateBlockData(StateBlockData): + CONFIG = ConfigBlock(implicit=True) + + def build(self): + super(NoPressureStateBlockData, self).build() + + self.flow_vol = Var(initialize=20) + self.flow_mol_phase_comp = Var( + self._params.phase_list, self._params.component_list, initialize=2 + ) + self.temperature = Var(initialize=300) + self.test_var = Var(initialize=1) + + def get_material_flow_terms(b, p, j): + return b.test_var + + def get_enthalpy_flow_terms(b, p): + return b.test_var + + def default_material_balance_type(self): + return MaterialBalanceType.componentPhase + + def default_energy_balance_type(self): + return EnergyBalanceType.enthalpyTotal + + +class TestMixer_NoPressure(object): + @declare_process_block_class("MixerFrame") + class MixerFrameData(MixerData): + def build(self): + super(MixerData, self).build() + + @pytest.mark.build + @pytest.mark.unit + def test_pressure_minimization_unsupported(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = NoPressureTestBlock() + + with pytest.raises( + PropertyNotSupportedError, + match="fs.scaler The property package supplied for this unit " + "does not appear to support pressure, which is required for " + "momentum mixing. Please set momentum_mixing_type to " + "MomentumMixingType.none or provide a property package which " + "supports pressure.", + ): + m.fs.scaler = Mixer( + property_package=m.fs.pp, + momentum_mixing_type=MomentumMixingType.minimize, + construct_ports=False, + ) + + @pytest.mark.build + @pytest.mark.unit + def test_pressure_equal_unsupported(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = NoPressureTestBlock() + + with pytest.raises( + PropertyNotSupportedError, + match="fs.scaler The property package supplied for this unit " + "does not appear to support pressure, which is required for " + "momentum mixing. Please set momentum_mixing_type to " + "MomentumMixingType.none or provide a property package which " + "supports pressure.", + ): + m.fs.scaler = Mixer( + property_package=m.fs.pp, + momentum_mixing_type=MomentumMixingType.equality, + construct_ports=False, + ) + + @pytest.mark.build + @pytest.mark.unit + def test_pressure_equal_and_min_unsupported(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = NoPressureTestBlock() + + with pytest.raises( + PropertyNotSupportedError, + match="fs.scaler The property package supplied for this unit " + "does not appear to support pressure, which is required for " + "momentum mixing. Please set momentum_mixing_type to " + "MomentumMixingType.none or provide a property package which " + "supports pressure.", + ): + m.fs.scaler = Mixer( + property_package=m.fs.pp, + momentum_mixing_type=MomentumMixingType.minimize_and_equality, + construct_ports=False, + ) + + @pytest.mark.build + @pytest.mark.unit + def test_pressure_none_unsupported(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = NoPressureTestBlock() + + m.fs.scaler = Mixer( + property_package=m.fs.pp, + momentum_mixing_type=MomentumMixingType.none, + construct_ports=False, + ) + + +# ----------------------------------------------------------------------------- +@pytest.mark.iapws +@pytest.mark.skipif(not iapws95.iapws95_available(), reason="IAPWS not available") +class TestIAPWS(object): + @pytest.fixture(scope="class") + def iapws(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = iapws95.Iapws95ParameterBlock() + + m.fs.unit = Mixer( + property_package=m.fs.properties, + material_balance_type=MaterialBalanceType.componentTotal, + momentum_mixing_type=MomentumMixingType.equality, + ) + + m.fs.unit.inlet_1.flow_mol[0].fix(100) + m.fs.unit.inlet_1.enth_mol[0].fix(5500) + m.fs.unit.inlet_1.pressure[0].fix(101325) + + m.fs.unit.inlet_2.flow_mol[0].fix(100) + m.fs.unit.inlet_2.enth_mol[0].fix(5000) + m.fs.unit.inlet_2.pressure[0].value = 1e5 + + return m + + @pytest.mark.build + @pytest.mark.unit + def test_build(self, iapws): + assert len(iapws.fs.unit.inlet_1.vars) == 3 + assert hasattr(iapws.fs.unit.inlet_1, "flow_mol") + assert hasattr(iapws.fs.unit.inlet_1, "enth_mol") + assert hasattr(iapws.fs.unit.inlet_1, "pressure") + + assert len(iapws.fs.unit.inlet_2.vars) == 3 + assert hasattr(iapws.fs.unit.inlet_2, "flow_mol") + assert hasattr(iapws.fs.unit.inlet_2, "enth_mol") + assert hasattr(iapws.fs.unit.inlet_2, "pressure") + + assert hasattr(iapws.fs.unit, "outlet") + assert len(iapws.fs.unit.outlet.vars) == 3 + assert hasattr(iapws.fs.unit.outlet, "flow_mol") + assert hasattr(iapws.fs.unit.outlet, "enth_mol") + assert hasattr(iapws.fs.unit.outlet, "pressure") + + @pytest.mark.component + def test_units(self, iapws): + assert_units_consistent(iapws) + + @pytest.mark.unit + def test_dof(self, iapws): + assert degrees_of_freedom(iapws) == 0 + + @pytest.mark.ui + @pytest.mark.unit + def test_get_stream_table_contents(self, iapws): + stable = iapws.fs.unit._get_stream_table_contents() + + expected = pandas.DataFrame.from_dict( + { + "Units": { + "Molar Flow": getattr(pyunits.pint_registry, "mole/second"), + "Mass Flow": getattr(pyunits.pint_registry, "kg/second"), + "T": getattr(pyunits.pint_registry, "K"), + "P": getattr(pyunits.pint_registry, "Pa"), + "Vapor Fraction": getattr(pyunits.pint_registry, "dimensionless"), + "Molar Enthalpy": getattr(pyunits.pint_registry, "J/mole"), + }, + "inlet_1": { + "Molar Flow": 100, + "Mass Flow": 1.8015, + "T": 346.05, + "P": 101325, + "Vapor Fraction": 0, + "Molar Enthalpy": 5500.0, + }, + "inlet_2": { + "Molar Flow": 100, + "Mass Flow": 1.8015, + "T": 339.43, + "P": 1e5, + "Vapor Fraction": 0, + "Molar Enthalpy": 5000, + }, + "Outlet": { + "Molar Flow": 1, + "Mass Flow": 1.8015e-2, + "T": 270.4877112932641, + "P": 11032305.8275, + "Vapor Fraction": 0, + "Molar Enthalpy": 0.01102138712926277, + }, + } + ) + + pandas.testing.assert_frame_equal(stable, expected, rtol=1e-4, atol=1e-4) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_initialize(self, iapws): + initialization_tester(iapws) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solve(self, iapws): + results = solver.solve(iapws) + + # Check for optimal solution + assert check_optimal_termination(results) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solution(self, iapws): + assert pytest.approx(200, abs=1e-5) == value(iapws.fs.unit.outlet.flow_mol[0]) + + assert pytest.approx(5250, abs=1e0) == value(iapws.fs.unit.outlet.enth_mol[0]) + + assert pytest.approx(101325, abs=1e2) == value(iapws.fs.unit.outlet.pressure[0]) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_conservation(self, iapws): + assert ( + abs( + value( + iapws.fs.unit.inlet_1.flow_mol[0] + + iapws.fs.unit.inlet_2.flow_mol[0] + - iapws.fs.unit.outlet.flow_mol[0] + ) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + iapws.fs.unit.inlet_1.flow_mol[0] + * iapws.fs.unit.inlet_1.enth_mol[0] + + iapws.fs.unit.inlet_2.flow_mol[0] + * iapws.fs.unit.inlet_2.enth_mol[0] + - iapws.fs.unit.outlet.flow_mol[0] + * iapws.fs.unit.outlet.enth_mol[0] + ) + ) + <= 1e-6 + ) + + +# ----------------------------------------------------------------------------- +class TestSaponification(object): + @pytest.fixture(scope="class") + def sapon(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = SaponificationParameterBlock() + + m.fs.unit = Mixer(property_package=m.fs.properties) + + m.fs.unit.inlet_1.flow_vol[0].fix(1e-3) + m.fs.unit.inlet_1.temperature[0].fix(320) + m.fs.unit.inlet_1.pressure[0].fix(101325) + m.fs.unit.inlet_1.conc_mol_comp[0, "H2O"].fix(55388.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "NaOH"].fix(100.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "EthylAcetate"].fix(100.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "SodiumAcetate"].fix(0.0) + m.fs.unit.inlet_1.conc_mol_comp[0, "Ethanol"].fix(0.0) + + m.fs.unit.inlet_2.flow_vol[0].fix(1e-3) + m.fs.unit.inlet_2.temperature[0].fix(300) + m.fs.unit.inlet_2.pressure[0].fix(101325) + m.fs.unit.inlet_2.conc_mol_comp[0, "H2O"].fix(55388.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "NaOH"].fix(100.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "EthylAcetate"].fix(100.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "SodiumAcetate"].fix(0.0) + m.fs.unit.inlet_2.conc_mol_comp[0, "Ethanol"].fix(0.0) + + return m + + @pytest.mark.build + @pytest.mark.unit + def test_build(self, sapon): + assert len(sapon.fs.unit.inlet_1.vars) == 4 + assert hasattr(sapon.fs.unit.inlet_1, "flow_vol") + assert hasattr(sapon.fs.unit.inlet_1, "conc_mol_comp") + assert hasattr(sapon.fs.unit.inlet_1, "temperature") + assert hasattr(sapon.fs.unit.inlet_1, "pressure") + + assert len(sapon.fs.unit.inlet_2.vars) == 4 + assert hasattr(sapon.fs.unit.inlet_2, "flow_vol") + assert hasattr(sapon.fs.unit.inlet_2, "conc_mol_comp") + assert hasattr(sapon.fs.unit.inlet_2, "temperature") + assert hasattr(sapon.fs.unit.inlet_2, "pressure") + + assert len(sapon.fs.unit.outlet.vars) == 4 + assert hasattr(sapon.fs.unit.outlet, "flow_vol") + assert hasattr(sapon.fs.unit.outlet, "conc_mol_comp") + assert hasattr(sapon.fs.unit.outlet, "temperature") + assert hasattr(sapon.fs.unit.outlet, "pressure") + + assert number_variables(sapon) == 26 + assert number_total_constraints(sapon) == 10 + assert number_unused_variables(sapon) == 0 + + @pytest.mark.component + def test_units(self, sapon): + assert_units_consistent(sapon) + + @pytest.mark.unit + def test_dof(self, sapon): + assert degrees_of_freedom(sapon) == 0 + + @pytest.mark.ui + @pytest.mark.unit + def test_get_stream_table_contents(self, sapon): + stable = sapon.fs.unit._get_stream_table_contents() + + expected = pandas.DataFrame.from_dict( + { + "Units": { + "Volumetric Flowrate": getattr( + pyunits.pint_registry, "m**3/second" + ), + "Molar Concentration H2O": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Molar Concentration NaOH": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Molar Concentration EthylAcetate": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Molar Concentration SodiumAcetate": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Molar Concentration Ethanol": getattr( + pyunits.pint_registry, "mole/m**3" + ), + "Temperature": getattr(pyunits.pint_registry, "K"), + "Pressure": getattr(pyunits.pint_registry, "Pa"), + }, + "inlet_1": { + "Volumetric Flowrate": 1e-3, + "Molar Concentration H2O": 55388, + "Molar Concentration NaOH": 100.00, + "Molar Concentration EthylAcetate": 100.00, + "Molar Concentration SodiumAcetate": 0, + "Molar Concentration Ethanol": 0, + "Temperature": 320, + "Pressure": 1.0132e05, + }, + "inlet_2": { + "Volumetric Flowrate": 1e-3, + "Molar Concentration H2O": 55388, + "Molar Concentration NaOH": 100.00, + "Molar Concentration EthylAcetate": 100.00, + "Molar Concentration SodiumAcetate": 0, + "Molar Concentration Ethanol": 0, + "Temperature": 300, + "Pressure": 1.0132e05, + }, + "Outlet": { + "Volumetric Flowrate": 1.00, + "Molar Concentration H2O": 100.00, + "Molar Concentration NaOH": 100.0, + "Molar Concentration EthylAcetate": 100.00, + "Molar Concentration SodiumAcetate": 100.00, + "Molar Concentration Ethanol": 100.00, + "Temperature": 298.15, + "Pressure": 1.0132e05, + }, + } + ) + + pandas.testing.assert_frame_equal(stable, expected, rtol=1e-4, atol=1e-4) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_initialize(self, sapon): + initialization_tester(sapon) + + @pytest.mark.component + def test_scaling(self, sapon): + sapon.fs.unit.calculate_scaling_factors() + + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solve(self, sapon): + results = solver.solve(sapon) + + # Check for optimal solution + assert check_optimal_termination(results) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_solution(self, sapon): + assert pytest.approx(2e-3, abs=1e-6) == value(sapon.fs.unit.outlet.flow_vol[0]) + + assert pytest.approx(55388.0, abs=1e0) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "H2O"] + ) + assert pytest.approx(100.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "NaOH"] + ) + assert pytest.approx(100.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "EthylAcetate"] + ) + assert pytest.approx(0.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "SodiumAcetate"] + ) + assert pytest.approx(0.0, abs=1e-3) == value( + sapon.fs.unit.outlet.conc_mol_comp[0, "Ethanol"] + ) + + assert pytest.approx(310.0, abs=1e-1) == value( + sapon.fs.unit.outlet.temperature[0] + ) + + assert pytest.approx(101325, abs=1e2) == value(sapon.fs.unit.outlet.pressure[0]) + + @pytest.mark.solver + @pytest.mark.skipif(solver is None, reason="Solver not available") + @pytest.mark.component + def test_conservation(self, sapon): + assert ( + abs( + value( + sapon.fs.unit.inlet_1.flow_vol[0] + + sapon.fs.unit.inlet_2.flow_vol[0] + - sapon.fs.unit.outlet.flow_vol[0] + ) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + sapon.fs.unit.inlet_1.flow_vol[0] + * sapon.fs.properties.dens_mol + * sapon.fs.properties.cp_mol + * ( + sapon.fs.unit.inlet_1.temperature[0] + - sapon.fs.properties.temperature_ref + ) + + sapon.fs.unit.inlet_2.flow_vol[0] + * sapon.fs.properties.dens_mol + * sapon.fs.properties.cp_mol + * ( + sapon.fs.unit.inlet_2.temperature[0] + - sapon.fs.properties.temperature_ref + ) + - sapon.fs.unit.outlet.flow_vol[0] + * sapon.fs.properties.dens_mol + * sapon.fs.properties.cp_mol + * ( + sapon.fs.unit.outlet.temperature[0] + - sapon.fs.properties.temperature_ref + ) + ) + ) + <= 1e-3 + ) + + +@pytest.mark.skipif(not cubic_roots_available(), reason="Cubic functions not available") +@pytest.mark.component +def test_construction_component_not_in_phase(): + m = ConcreteModel() + m.fs = FlowsheetBlock() + m.fs.prop_params = GenericParameterBlock(**get_prop(["H2O", "H2"], ["Liq", "Vap"])) + m.fs.inject1 = Mixer( + property_package=m.fs.prop_params, + inlet_list=["in1", "in2"], + momentum_mixing_type=MomentumMixingType.none, + ) + iscale.calculate_scaling_factors(m) + + +@pytest.mark.unit +def test_initialization_error(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.pp = PhysicalParameterTestBlock() + + m.fs.scaler = Mixer(property_package=m.fs.pp) + + m.fs.scaler.inlet_1_state[0].material_flow_mol.fix(10) + m.fs.scaler.inlet_2_state[0].material_flow_mol.fix(10) + m.fs.scaler.mixed_state[0].material_flow_mol.fix(100) + + with pytest.raises(InitializationError): + m.fs.scaler.initialize() diff --git a/idaes/models_extra/column_models/MEAsolvent_column.py b/idaes/models_extra/column_models/MEAsolvent_column.py index ff84e57863..b3051530ba 100644 --- a/idaes/models_extra/column_models/MEAsolvent_column.py +++ b/idaes/models_extra/column_models/MEAsolvent_column.py @@ -1,14 +1,14 @@ -################################################################################# +################################################################################ # 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). +# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021 +# 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. # -# Copyright (c) 2018-2023 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. +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and +# license information. ################################################################################# """ Packed Solvent Column Model for MEA systems @@ -23,17 +23,31 @@ import numpy as np # Import Pyomo libraries +import pyomo.opt from pyomo.environ import ( + ConcreteModel, value, Var, NonNegativeReals, + Param, + TransformationFactory, Constraint, Expression, + Objective, + SolverStatus, + TerminationCondition, check_optimal_termination, + assert_optimal_termination, exp, + log, units as pyunits, + Set, + Reference, ) +from pyomo.common.collections import ComponentSet, ComponentMap + from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.common.config import ConfigValue, Bool # Import IDAES Libraries from idaes.core.util.constants import Constants as CONST @@ -45,10 +59,26 @@ from idaes.core.solvers.get_solver import get_solver import idaes.logger as idaeslog +from idaes.core.solvers import use_idaes_solver_configuration_defaults import idaes.core.util.scaling as iscale - - -__author__ = "Paul Akula, John Eslick, Anuja Deshpande, Andrew Lee" +from pyomo.util.subsystems import ( + create_subsystem_block, +) +from idaes.core.solvers.petsc import ( + _sub_problem_scaling_suffix, +) +from idaes.core.initialization import BlockTriangularizationInitializer +from idaes.core.util.initialization import _fix_vars, _restore_fixedness +from idaes.models_extra.column_models.enhancement_factor_model_pseudo_second_order import ( + make_enhancement_factor_model, + initialize_enhancement_factor_model +) +# from idaes.models_extra.column_models.enhancement_factor_model_third_order import ( +# make_enhancement_factor_model, +# initialize_enhancement_factor_model +# ) +from idaes.core.surrogate.surrogate_block import SurrogateBlock +__author__ = "Paul Akula, John Eslick, Anuja Deshpande, Andrew Lee, Douglas Allan" @declare_process_block_class("MEAColumn") @@ -59,14 +89,44 @@ class MEAColumnData(PackedColumnData): CONFIG = PackedColumnData.CONFIG() + CONFIG.declare( + "corrected_hx_coeff_eqn", + ConfigValue( + default=False, + domain=Bool, + description="Use heat transfer equation appropriate for large material fluxes", + doc="""Boolean flag indicating whether to use correction factor for heat transfer + equation appropriate for when 'large' material transfer fluxes are in column. + Warning: using the correction factor when fluxes are 'small' can result in badly + conditioned problem.""", + ), + ) + CONFIG.declare( + "surrogate_enhancement_factor_model", + ConfigValue( + default=None, + description="Surrogate object for enhancement factor model. Use first-principles model if None.", + doc="""Surrogate object to use for enhancement factor model. Must take six inputs: CO2_loading and + H2O_loading of the liquid stream, the logarithm of the partial pressure of CO2 in the vapor stream, + the temperature of the liquid stream, the liquid mass transfer coefficient, and the vapor mass transfer + coefficient, in that order. One output, the logarithm of the enhancement factor. + """, + ), + ) + def liquid_phase_mass_transfer_model(self): """ Enhancement factor based liquid phase mass transfer model. """ + time = self.flowsheet().time vap_comp = self.config.vapor_phase.property_package.component_list liq_comp = self.config.liquid_phase.property_package.component_list equilibrium_comp = vap_comp & liq_comp solute_comp_list = ["CO2"] + log_diffus_liq_comp_list = self.log_diffus_liq_comp_list = [ + "CO2", + "MEA", + ] # Can add ions if we want them lunits = ( self.config.liquid_phase.property_package.get_metadata().get_derived_units @@ -77,38 +137,47 @@ def liquid_phase_mass_transfer_model(self): self.flowsheet().time, self.liquid_phase.length_domain, solute_comp_list, + bounds=(0, None), units=lunits("length") / lunits("time"), doc="Liquid phase mass transfer coefficient", ) - self.enhancement_factor = Var( + self.log_enhancement_factor = Var( self.flowsheet().time, self.liquid_phase.length_domain, units=pyunits.dimensionless, - initialize=160, + bounds=(None, 100), # 100 is about where we start getting AMPL evaluation errors due to overflow + initialize=5, + doc="Natural logarithm of the enhancement factor", + ) + + @self.Expression( + self.flowsheet().time, + self.liquid_phase.length_domain, doc="Enhancement factor", ) + def enhancement_factor(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return exp(blk.log_enhancement_factor[t, x]) - # Intermediate term - def rule_phi(blk, t, x, j): - if x == self.vapor_phase.length_domain.first(): + @self.Expression( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Intermediate for calculating CO2 mass transfer driving force", + ) + def psi(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): return Expression.Skip else: zb = self.liquid_phase.length_domain.prev(x) return ( blk.enhancement_factor[t, zb] - * blk.mass_transfer_coeff_liq[t, zb, j] - / blk.mass_transfer_coeff_vap[t, x, j] + * blk.mass_transfer_coeff_liq[t, zb, "CO2"] + / blk.mass_transfer_coeff_vap[t, x, "CO2"] ) - self.phi = Expression( - self.flowsheet().time, - self.vapor_phase.length_domain, - solute_comp_list, - rule=rule_phi, - doc="Equilibrium partial pressure intermediate term", - ) - @self.Constraint( self.flowsheet().time, self.vapor_phase.length_domain, @@ -116,1122 +185,1420 @@ def rule_phi(blk, t, x, j): doc="Equilibrium partial pressure at interface", ) def pressure_at_interface(blk, t, x, j): - if x == self.vapor_phase.length_domain.first(): + if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: zb = self.liquid_phase.length_domain.prev(x) lprops = blk.liquid_phase.properties[t, zb] henrycomp = lprops.params.get_component(j).config.henry_component + Pressure = pyunits.convert( + blk.vapor_phase.properties[t, x].pressure, + to_units=lunits("pressure"), + ) if henrycomp is not None and "Liq" in henrycomp: - return blk.pressure_equil[t, x, j] == ( - ( - blk.vapor_phase.properties[t, x].mole_frac_comp[j] - * pyunits.convert( - blk.vapor_phase.properties[t, x].pressure, - to_units=lunits("pressure"), - ) - + blk.phi[t, x, j] - * lprops.conc_mol_phase_comp_true["Liq", j] - ) + return blk.mass_transfer_driving_force[t, x, j] == ( + blk.psi[t, x] / ( - 1 - + blk.phi[t, x, j] - / blk.liquid_phase.properties[t, zb].henry["Liq", j] + blk.psi[t, x] + + blk.liquid_phase.properties[t, zb].henry["Liq", j] + ) + * ( + blk.vapor_phase.properties[t, x].mole_frac_comp[j] + * Pressure + - lprops.conc_mol_phase_comp_true["Liq", j] + * blk.liquid_phase.properties[t, zb].henry["Liq", j] ) ) else: - return blk.pressure_equil[t, x, j] == ( - lprops.mole_frac_phase_comp_true["Liq", j] + return blk.mass_transfer_driving_force[t, x, j] == ( + blk.vapor_phase.properties[t, x].mole_frac_comp[j] * Pressure + - lprops.mole_frac_phase_comp_true["Liq", j] * lprops.pressure_sat_comp[j] ) def build(self): - super().build() # --------------------------------------------------------------------- # Unit level sets vap_comp = self.config.vapor_phase.property_package.component_list liq_comp = self.config.liquid_phase.property_package.component_list - equilibrium_comp = vap_comp & liq_comp - solute_comp_list = ["CO2"] + self.equilibirum_comp = equilibrium_comp = vap_comp & liq_comp + solute_comp_list = self.solute_comp_list = Set( + initialize=["CO2"], + ordered=True, + doc="Set of solutes", + ) lunits = ( self.config.liquid_phase.property_package.get_metadata().get_derived_units ) - # Velocity calculations - self.velocity_vap = Var( + # Map from log variable of property to equation that defines it. + self.log_property_var_eqn_map = ComponentMap() + # Map from log variable of parameter (fixed var) to equation that defines it. + self.log_parameter_var_eqn_map = ComponentMap() + + # Log variables from property package + self.log_dens_mass_liq = Var( self.flowsheet().time, - self.vapor_phase.length_domain, - domain=NonNegativeReals, - initialize=2, - units=pyunits.m / pyunits.s, - doc="Vapor superficial velocity", + self.liquid_phase.length_domain, + initialize=7, + units=pyunits.dimensionless, + doc="Logarithm of liquid mass density", ) - self.velocity_liq = Var( + @self.Constraint( self.flowsheet().time, self.liquid_phase.length_domain, - domain=NonNegativeReals, - initialize=0.01, - units=pyunits.m / pyunits.s, - doc="Liquid superficial velocity", + doc="Constraint defining log density variable", ) - - def eq_velocity_vap(blk, t, x): + def log_dens_mass_liq_eqn(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip return ( - blk.velocity_vap[t, x] - * blk.area_column - * blk.vapor_phase.properties[t, x].dens_mol - == blk.vapor_phase.properties[t, x].flow_mol + exp(blk.log_dens_mass_liq[t, x]) * lunits("density_mass") + == blk.liquid_phase.properties[t, x].dens_mass_phase["Liq"] ) - self.velocity_vap_eqn = Constraint( + self.log_property_var_eqn_map[ + self.log_dens_mass_liq + ] = self.log_dens_mass_liq_eqn + + self.log_surf_tens_liq = Var( self.flowsheet().time, - self.vapor_phase.length_domain, - rule=eq_velocity_vap, - doc="Vapor superficial velocity", + self.liquid_phase.length_domain, + initialize=-3, + units=pyunits.dimensionless, + doc="Logarithm of liquid surface tension", ) - def eq_velocity_liq(blk, t, x): + @self.Constraint( + self.flowsheet().time, + self.liquid_phase.length_domain, + doc="Constraint defining log surface tension variable", + ) + def log_surf_tens_liq_eqn(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip return ( - blk.velocity_liq[t, x] - * blk.area_column - * blk.liquid_phase.properties[t, x].dens_mol - == blk.liquid_phase.properties[t, x].flow_mol + exp(blk.log_surf_tens_liq[t, x]) * lunits("surface_tension") + == blk.liquid_phase.properties[t, x].surf_tens_phase["Liq"] ) - self.velocity_liq_eqn = Constraint( + self.log_property_var_eqn_map[ + self.log_surf_tens_liq + ] = self.log_surf_tens_liq_eqn + + self.log_visc_d_liq = Var( self.flowsheet().time, self.liquid_phase.length_domain, - rule=eq_velocity_liq, - doc="Liquid superficial velocity", + bounds=(None, 100), + initialize=1, + units=pyunits.dimensionless, + doc="Logarithm of liquid viscosity", ) - # --------------------------------------------------------------------- - # Vapor-liquid heat transfer coeff modified by Ackmann factor - self.heat_transfer_coeff_base = Var( + @self.Constraint( self.flowsheet().time, - self.vapor_phase.length_domain, - initialize=100, - units=lunits("power") / lunits("temperature") / lunits("length") ** 2, - doc="Uncorrected vapor-liquid heat transfer coefficient", + self.liquid_phase.length_domain, + doc="Defines log variable for liquid viscosity", ) - - def rule_heat_transfer_coeff_Ack(blk, t, x): - if x == self.vapor_phase.length_domain.first(): + def log_visc_d_liq_eqn(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): return Constraint.Skip else: - Ackmann_factor = sum( - blk.vapor_phase.properties[t, x].cp_mol_phase_comp["Vap", j] - * blk.interphase_mass_transfer[t, x, j] - for j in equilibrium_comp - ) - return blk.heat_transfer_coeff[t, x] == Ackmann_factor / ( - 1 - - exp( - -Ackmann_factor - / ( - blk.heat_transfer_coeff_base[t, x] - * blk.area_interfacial[t, x] - * blk.area_column - ) - ) + return ( + exp(blk.log_visc_d_liq[t, x]) * lunits("dynamic_viscosity") + == blk.liquid_phase.properties[t, x].visc_d_phase["Liq"] ) - self.heat_transfer_coeff_corr = Constraint( - self.flowsheet().time, - self.vapor_phase.length_domain, - rule=rule_heat_transfer_coeff_Ack, - doc="Vap-Liq heat transfer correction by Ackmann factor", - ) - - # Interfacial Area model ([m2/m3]): - # Reference: Tsai correlation,regressed by Chinen et al. 2018 - - self.area_interfacial_parA = Var( - initialize=0.6486, units=None, doc="Interfacial area parameter A" - ) + self.log_property_var_eqn_map[self.log_visc_d_liq] = self.log_visc_d_liq_eqn - self.area_interfacial_parB = Var( - initialize=0.12, + self.log_diffus_liq_comp = Var( + self.flowsheet().time, + self.liquid_phase.length_domain, + self.log_diffus_liq_comp_list, + bounds=(None, 100), + initialize=1, units=pyunits.dimensionless, - doc="Interfacial area parameter B", + doc="""Logarithm of the liquid phase diffusivity of a species""", ) - def rule_interfacial_area(blk, t, x): - if x == self.vapor_phase.length_domain.first(): + @self.Constraint( + self.flowsheet().time, + self.liquid_phase.length_domain, + self.log_diffus_liq_comp_list, + doc="Defines log variable for liquid phase diffusivity", + ) + def log_diffus_liq_comp_eqn(blk, t, x, j): + if x == blk.liquid_phase.length_domain.last(): return Constraint.Skip else: - return ( - blk.area_interfacial[t, x] - == blk.packing_specific_area - * blk.area_interfacial_parA - * ( - ( - pyunits.newton - * (pyunits.m ** (2 / 3)) - * (pyunits.s ** (4 / 3)) - / pyunits.kilogram - ) - ** blk.area_interfacial_parB - ) - * ( - ( - blk.liquid_phase.properties[t, x].mw - / blk.liquid_phase.properties[t, x].vol_mol_phase["Liq"] - / blk.liquid_phase.properties[t, x].surf_tens_phase["Liq"] - ) - * (blk.velocity_liq[t, x]) ** (4.0 / 3.0) - ) - ** blk.area_interfacial_parB - ) + return exp(blk.log_diffus_liq_comp[t, x, j]) * lunits( + "diffusivity" + ) == (blk.liquid_phase.properties[t, x].diffus_phase_comp["Liq", j]) + + self.log_property_var_eqn_map[ + self.log_diffus_liq_comp + ] = self.log_diffus_liq_comp_eqn - self.area_interfacial_constraint = Constraint( + self.log_dens_mass_vap = Var( self.flowsheet().time, self.vapor_phase.length_domain, - rule=rule_interfacial_area, - doc="Specific interfacial area", + bounds=(None, 100), + initialize=1, + units=pyunits.dimensionless, + doc="Logarithm of the vapor phase mass density", ) - # Liquid holdup model - # Reference: Tsai correlation,regressed by Chinen et al. 2018 + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def log_dens_mass_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + return exp(blk.log_dens_mass_vap[t, x]) * lunits( + "density_mass" + ) == pyunits.convert( + blk.vapor_phase.properties[t, x].dens_mass_phase["Vap"], + to_units=lunits("density_mass"), + ) - self.holdup_parA = Var(initialize=24.23, units=None, doc="Holdup parameter A") + self.log_property_var_eqn_map[ + self.log_dens_mass_vap + ] = self.log_dens_mass_vap_eqn - self.holdup_parB = Var( - initialize=0.6471, units=pyunits.dimensionless, doc="Holdup parameter B" + self.log_visc_d_vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + initialize=-10, + units=pyunits.dimensionless, + doc="Logarithm of vapor viscosity", ) - def rule_holdup_liq(blk, t, x): - if x == self.liquid_phase.length_domain.last(): + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Defines log variable for vapor viscosity", + ) + def log_visc_d_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: - return ( - blk.holdup_liq[t, x] - == blk.holdup_parA - * ( - ( - (pyunits.pascal * (pyunits.s ** (10 / 3))) - / (pyunits.kilogram * (pyunits.m ** (2 / 3))) - ) - ** blk.holdup_parB - ) - * ( - blk.velocity_liq[t, x] - * ( - blk.liquid_phase.properties[t, x].visc_d_phase["Liq"] - / ( - blk.liquid_phase.properties[t, x].mw - / blk.liquid_phase.properties[t, x].vol_mol_phase["Liq"] - ) - ) - ** (1 / 3) - ) - ** blk.holdup_parB + return exp(blk.log_visc_d_vap[t, x]) * lunits( + "dynamic_viscosity" + ) == pyunits.convert( + blk.vapor_phase.properties[t, x].visc_d_phase["Vap"], + to_units=lunits("dynamic_viscosity"), ) - self.holdup_liq_constraint = Constraint( + self.log_property_var_eqn_map[self.log_visc_d_vap] = self.log_visc_d_vap_eqn + + self.log_diffus_vap_comp = Var( self.flowsheet().time, - self.liquid_phase.length_domain, - rule=rule_holdup_liq, - doc="Volumetric liquid holdup [-]", + self.vapor_phase.length_domain, + equilibrium_comp, + initialize=1, + units=pyunits.dimensionless, + doc="""Logarithm of the vapor mass transfer coeff""", ) - # Mass transfer coefficients of diffusing components in vapor phase [mol/m2.s.Pa] - # Mass transfer coefficients, Billet and Schultes (1999) correlation, - # where parameters are regressed by Chinen et al. (2018). - - self.Cv_ref = Var( - initialize=0.357, - doc="Vapor packing specific constant in Billet and Schultes " - "volumetric mass transfer coefficient correlation", + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + equilibrium_comp, + doc="""Defines log variable for the vapor species diffusivity""", ) - - def rule_mass_transfer_coeff_vap(blk, t, x, j): - if x == self.vapor_phase.length_domain.first(): + def log_diffus_vap_comp_eqn(blk, t, x, j): + if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: - return blk.mass_transfer_coeff_vap[t, x, j] == ( - 1 - / ( - pyunits.convert( - CONST.gas_constant, - pyunits.kilogram - * pyunits.m**2 - / (pyunits.s**2 * pyunits.mol * pyunits.K), - ) - * blk.vapor_phase.properties[t, x].temperature - ) - ) * (blk.Cv_ref / (blk.holdup_vap[t, x]) ** 0.5) * ( - (blk.packing_specific_area / blk.hydraulic_diameter) ** 0.5 - ) * ( - (blk.vapor_phase.properties[t, x].diffus_phase_comp["Vap", j]) - ** (2 / 3) - ) * ( - ( - blk.vapor_phase.properties[t, x].visc_d_phase["Vap"] - / ( - blk.vapor_phase.properties[t, x].mw - / blk.vapor_phase.properties[t, x].vol_mol_phase["Vap"] - ) - ) - ** (1 / 3) - ) * ( - ( - ( - blk.velocity_vap[t, x] - * ( - blk.vapor_phase.properties[t, x].mw - / blk.vapor_phase.properties[t, x].vol_mol_phase["Vap"] - ) - ) - / ( - blk.packing_specific_area - * blk.vapor_phase.properties[t, x].visc_d_phase["Vap"] - ) - ) - ** (3 / 4) + diffus_vap_comp = pyunits.convert( + blk.vapor_phase.properties[t, x].diffus_phase_comp["Vap", j], + to_units=lunits("diffusivity"), + ) + return ( + exp(blk.log_diffus_vap_comp[t, x, j]) * lunits("diffusivity") + == diffus_vap_comp ) - self.mass_transfer_coeff_vap_constraint = Constraint( + self.log_property_var_eqn_map[ + self.log_diffus_vap_comp + ] = self.log_diffus_vap_comp_eqn + + self.log_pressure_vap = Var( self.flowsheet().time, self.vapor_phase.length_domain, - equilibrium_comp, - rule=rule_mass_transfer_coeff_vap, - doc=" Vapor phase mass transfer coefficient", + initialize=11.5, + units=pyunits.dimensionless, ) - # Mass transfer coefficients of diffusing components in liquid phase [m/s] - # Mass transfer coefficients, Billet and Schultes (1999) correlation, - # where parameters are regressed by Chinen et al. (2018). + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def log_pressure_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + pressure = pyunits.convert( + blk.vapor_phase.properties[t, x].pressure, + to_units=lunits("pressure"), + ) + return exp(blk.log_pressure_vap[t, x]) * lunits("pressure") == pressure - self.Cl_ref = Var( - initialize=0.5, - doc="""Liquid packing specific constant in Billet and Schultes - volumetric mass transfer coefficient correlation""", + self.log_property_var_eqn_map[self.log_pressure_vap] = self.log_pressure_vap_eqn + + self.log_dens_mol_vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + initialize=0, + units=pyunits.dimensionless, ) - def rule_mass_transfer_coeff_liq(blk, t, x, j): - if x == self.liquid_phase.length_domain.last(): + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def log_dens_mol_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: + dens_mol_vap = pyunits.convert( + blk.vapor_phase.properties[t, x].dens_mol_phase["Vap"], + to_units=lunits("density_mole"), + ) return ( - blk.mass_transfer_coeff_liq[t, x, j] - == blk.Cl_ref - * (12 ** (1 / 6)) - * ( - blk.velocity_liq[t, x] - * blk.liquid_phase.properties[t, x].diffus_phase_comp["Liq", j] - / (blk.hydraulic_diameter * blk.holdup_liq[t, x]) - ) - ** 0.5 + exp(blk.log_dens_mol_vap[t, x]) * lunits("density_mole") + == dens_mol_vap ) - self.mass_transfer_coeff_liq_constraint = Constraint( + self.log_property_var_eqn_map[self.log_dens_mol_vap] = self.log_dens_mol_vap_eqn + + self.log_therm_cond_vap = Var( self.flowsheet().time, - self.liquid_phase.length_domain, - solute_comp_list, - rule=rule_mass_transfer_coeff_liq, - doc="Liquid phase mass transfer coefficient", + self.vapor_phase.length_domain, + initialize=-3.5, + units=pyunits.dimensionless, ) - # Fix mass transfer parameters - self.Cv_ref.fix() - self.Cl_ref.fix() - - # Fix interfacial area parameters - self.area_interfacial_parA.fix() - self.area_interfacial_parB.fix() - - # Fix liquid holdup parameters - self.holdup_parA.fix() - self.holdup_parB.fix() + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def log_therm_cond_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + therm_cond_vap = pyunits.convert( + blk.vapor_phase.properties[t, x].therm_cond_phase["Vap"], + to_units=lunits("thermal_conductivity"), + ) + return ( + exp(blk.log_therm_cond_vap[t, x]) * lunits("thermal_conductivity") + == therm_cond_vap + ) - # CO2 molar concentration at interface - def rule_co2_conc_interface(blk, t, x): + self.log_property_var_eqn_map[ + self.log_therm_cond_vap + ] = self.log_therm_cond_vap_eqn - if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip + self.log_cp_mol_vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + initialize=3.5, + units=pyunits.dimensionless, + ) + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def log_cp_mol_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip else: - zf = blk.liquid_phase.length_domain.next(x) - + cp_mol_vap = pyunits.convert( + blk.vapor_phase.properties[t, x].cp_mol_phase["Vap"], + to_units=lunits("heat_capacity_mole"), + ) return ( - blk.pressure_equil[t, zf, "CO2"] - / blk.liquid_phase.properties[t, x].henry["Liq", "CO2"] + exp(blk.log_cp_mol_vap[t, x]) * lunits("heat_capacity_mole") + == cp_mol_vap ) - self.co2_conc_interface = Expression( + self.log_property_var_eqn_map[self.log_cp_mol_vap] = self.log_cp_mol_vap_eqn + + # Velocity calculations + self.velocity_vap = Var( self.flowsheet().time, - self.liquid_phase.length_domain, - rule=rule_co2_conc_interface, - doc="""CO2 concentration at interface""", + self.vapor_phase.length_domain, + domain=NonNegativeReals, + initialize=2, + units=lunits("velocity"), + doc="Vapor superficial velocity", ) - # Enhancement factor model - # Reference: Jozsef Gaspar,Philip Loldrup Fosbol, (2015) - - self.conc_interface_MEA = Var( + self.velocity_liq = Var( self.flowsheet().time, self.liquid_phase.length_domain, - bounds=(0.5, 1), - initialize=1, - units=pyunits.dimensionless, - doc="""Dimensionless concentration of MEA - at interface """, + domain=NonNegativeReals, + initialize=0.01, + units=pyunits.m / pyunits.s, + doc="Liquid superficial velocity", ) - self.sqrt_conc_interface_MEA = Var( + @self.Constraint( self.flowsheet().time, - self.liquid_phase.length_domain, - bounds=(0.5, 1), - initialize=0.97, - units=pyunits.dimensionless, - doc="""Substitute for conc_interface_MEA""", + self.vapor_phase.length_domain, + doc="Vapor superficial velocity", ) + def velocity_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + return ( + blk.velocity_vap[t, x] + * blk.area_column + * blk.vapor_phase.properties[t, x].dens_mol + == blk.vapor_phase.properties[t, x].flow_mol + ) - def rule_rate_constant(blk, t, x): + @self.Constraint( + self.flowsheet().time, + self.liquid_phase.length_domain, + doc="Liquid superficial velocity", + ) + def velocity_liq_eqn(blk, t, x): if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip - else: - T = blk.liquid_phase.properties[t, x].temperature - C_MEA = blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEA" - ] - C_H2O = blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "H2O" - ] + return Constraint.Skip + return ( + blk.velocity_liq[t, x] + * blk.area_column + * blk.liquid_phase.properties[t, x].dens_mol + == blk.liquid_phase.properties[t, x].flow_mol + ) - # Reference: X.Luo et al., Chem. Eng. Sci. (2015) + # Log-form variables and constraints for vapor and liquid velocities + # used in reformulations of correlations + self.log_velocity_vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=0.69, + doc="Logarithm of vapor velocity", + ) + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Defines log_velocity_vap", + ) + def log_velocity_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: return ( - ( - 2.003e10 * exp(-4742 * pyunits.K / T) * C_MEA - + 4.147e6 * exp(-3110 * pyunits.K / T) * C_H2O - ) - * 1e-6 - * ((pyunits.m) ** 6 / (pyunits.mol**2 * pyunits.s)) + exp(blk.log_velocity_vap[t, x]) * lunits("velocity") + == blk.velocity_vap[t, x] ) - self.rate_constant = Expression( + self.log_velocity_liq = Var( self.flowsheet().time, self.liquid_phase.length_domain, - rule=rule_rate_constant, - doc="Second order rate constant [m3/(mol.s)]", + bounds=(None, 100), + initialize=-4.6, + doc="""Logarithm of liquid velocity""", ) - def rule_hatta_number(blk, t, x): - - if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip - else: - return ( - blk.rate_constant[t, x] - * ( - blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEA" - ] - ) - * blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ - "Liq", "CO2" - ] - ) ** 0.5 / blk.mass_transfer_coeff_liq[t, x, "CO2"] - - self.Hatta = Expression( + @self.Constraint( self.flowsheet().time, self.liquid_phase.length_domain, - rule=rule_hatta_number, - doc="""Hatta number""", + doc="Defines log_velocity_liq", ) - - def rule_conc_CO2_bulk(blk, t, x): + def log_velocity_liq_eqn(blk, t, x): if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip + return Constraint.Skip else: return ( - blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "CO2" - ] - / blk.co2_conc_interface[t, x] + exp(blk.log_velocity_liq[t, x]) * lunits("velocity") + == blk.velocity_liq[t, x] ) - self.conc_CO2_bulk = Expression( + self.log_property_var_eqn_map[self.log_velocity_liq] = self.log_velocity_liq_eqn + self.log_property_var_eqn_map[self.log_velocity_vap] = self.log_velocity_vap_eqn + + @self.Expression( self.flowsheet().time, - self.liquid_phase.length_domain, - rule=rule_conc_CO2_bulk, - doc="""Dimensionless concentration of CO2, - Driving force term where - Absorption implies conc_CO2_bulk < 1 and - Desorption impies conc_CO2_bulk > 1 """, + self.vapor_phase.length_domain, + doc="Partial pressure of CO2 Henry's law [Pa]", ) - - def rule_instantaneous_E(blk, t, x): - if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip + def PpCO2_equil_He(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return 1000 else: - return 1 + ( - blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ - "Liq", "MEA" - ] - * blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEA" - ] - ) / ( - 2 - * blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ + zb = self.liquid_phase.length_domain.prev(x) + return ( + blk.liquid_phase.properties[t, zb].conc_mol_phase_comp_true[ "Liq", "CO2" ] - * blk.co2_conc_interface[t, x] + * blk.liquid_phase.properties[t, zb].henry["Liq", "CO2"] ) - self.instant_E = Expression( + # ##################################################################### + # Interfacial Area model ([m2/m3]): + # Reference: Tsai correlation,regressed by Chinen et al. 2018 + + @self.Expression() + def wetted_perimeter(blk): + return blk.area_column * blk.packing_specific_area / blk.eps_ref + + # Original value with weird fractional power units was 0.6486 + # Corrected value was calculated by unlumping parameters to + # permit for unit handling + self.area_interfacial_parA = Var( + initialize=1.43914, + units=pyunits.dimensionless, + doc="Interfacial area parameter A", + ) + + self.area_interfacial_parB = Var( + initialize=0.12, + units=pyunits.dimensionless, + doc="Interfacial area parameter B", + ) + + self.log_area_interfacial = Var( self.flowsheet().time, - self.liquid_phase.length_domain, - rule=rule_instantaneous_E, - doc="Instantaneous Enhancement factor", + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=5, + doc="Logarithm of interfacial area", ) - def rule_conc_interface_MEACOO(blk, t, x): - if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Defines log variable for interfacial area", + ) + def log_area_interfacial_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip else: - return 1 + ( - blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ - "Liq", "MEA" - ] - * blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEA" - ] - ) * ( - 1 - - blk.sqrt_conc_interface_MEA[t, x] - * blk.sqrt_conc_interface_MEA[t, x] - ) / ( - 2 - * blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ - "Liq", "MEACOO_-" - ] - * blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEACOO_-" - ] + return ( + exp(blk.log_area_interfacial[t, x]) + * lunits("length") ** 2 + / lunits("length") ** 3 + == blk.area_interfacial[t, x] ) - self.conc_interface_MEACOO = Expression( - self.flowsheet().time, - self.liquid_phase.length_domain, - rule=rule_conc_interface_MEACOO, - doc="Dimensionless concentration of MEACOO-", + self.log_area_interfacial_parA = Var( + bounds=(None, 100), + initialize=0.3640457058674088, + doc="Logarithm of interfacial area model param A", ) - def rule_conc_interface_MEAH(blk, t, x): - if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip + @self.Constraint(doc="Defines log variable for interfacial area param A") + def log_area_interfacial_parA_eqn(blk): + return exp(blk.log_area_interfacial_parA) == blk.area_interfacial_parA + + self.log_parameter_var_eqn_map[ + self.log_area_interfacial_parA + ] = self.log_area_interfacial_parA_eqn + + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Defines specific interfacial area", + ) + def area_interfacial_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip else: - return 1 + ( - blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ - "Liq", "MEA" - ] - * blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEA" - ] - ) * ( - 1 - - blk.sqrt_conc_interface_MEA[t, x] - * blk.sqrt_conc_interface_MEA[t, x] - ) / ( - 2 - * blk.liquid_phase.properties[t, x].diffus_phase_comp_true[ - "Liq", "MEA_+" - ] - * blk.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ - "Liq", "MEA_+" - ] + x_liq = blk.liquid_phase.length_domain.prev(x) + log_packing_specific_area = log( + blk.packing_specific_area + / (lunits("length") ** 2 / lunits("length") ** 3) + ) + log_g = log( + pyunits.convert( + CONST.acceleration_gravity, to_units=lunits("acceleration") + ) + / lunits("acceleration") + ) + log_cross_sectional_area = log(blk.area_column / lunits("area")) + log_wetted_perimeter = log(blk.wetted_perimeter / lunits("length")) + return blk.log_area_interfacial[t, x] == ( + log_packing_specific_area + + blk.log_area_interfacial_parA + + blk.area_interfacial_parB + * ( + blk.log_dens_mass_liq[t, x_liq] + - blk.log_surf_tens_liq[t, x_liq] + + (1 / 3) * log_g + + (4 / 3) + * ( + blk.log_velocity_liq[t, x_liq] + + log_cross_sectional_area + - log_wetted_perimeter + ) + ) ) - self.conc_interface_MEAH = Expression( + # ##################################################################### + # Liquid holdup model + # Reference: Tsai correlation,regressed by Chinen et al. 2018 + + # Original, dimensional value is 24.23. Corrected by subtracting the contribution of alpha parameter + # to make this parameter dimensionless + self.holdup_parA = Var( + initialize=11.4474, units=pyunits.dimensionless, doc="Holdup parameter A" + ) + + self.holdup_parB = Var( + initialize=0.6471, units=pyunits.dimensionless, doc="Holdup parameter B" + ) + + # Value taken from CCSI model. Most of the factors lumped together into this factor can be calculated + # from information we have, but the channel size appears only in this equation for the set of correlations + # we are using. Technically that parameter exists in the SolventColumn model, but we have nothing on which + # to base its value. + # TODO Bring definition of these parameters into flowsheets + self.holdup_parAlpha = Param( + initialize=3.185966, + units=1 / (lunits("length") * (lunits("acceleration") ** (2 / 3))), + doc="Holdup parameter alpha", + mutable=True, + ) + + self.log_holdup_liq = Var( self.flowsheet().time, self.liquid_phase.length_domain, - rule=rule_conc_interface_MEAH, - doc="Dimensionless concentration of MEA+", + bounds=(None, 100), + initialize=-1, + units=pyunits.dimensionless, + doc="Logarithm of liquid holdup", ) - @self.Expression( + @self.Constraint( self.flowsheet().time, self.liquid_phase.length_domain, - doc="""Dimensionless concentration of CO2 - at equilibrium with the bulk """, + doc="Defines log variable for liquid holdup", ) - def conc_CO2_equil_bulk(blk, t, x): + def log_holdup_liq_eqn(blk, t, x): if x == blk.liquid_phase.length_domain.last(): - return Expression.Skip + return Constraint.Skip else: - return ( - blk.conc_CO2_bulk[t, x] - * blk.conc_interface_MEAH[t, x] - * blk.conc_interface_MEACOO[t, x] - / (blk.sqrt_conc_interface_MEA[t, x] ** 4) - ) + return exp(blk.log_holdup_liq[t, x]) == blk.holdup_liq[t, x] + + self.log_holdup_parA = Var( + bounds=(None, 100), + initialize=3.1875915348284343, + doc="Logarithm of liquid holdup model param A", + ) + + @self.Constraint(doc="Defines log variable for liquid holdup param A") + def log_holdup_parA_eqn(blk): + return exp(blk.log_holdup_parA) == blk.holdup_parA + + self.log_parameter_var_eqn_map[self.log_holdup_parA] = self.log_holdup_parA_eqn @self.Constraint( self.flowsheet().time, self.liquid_phase.length_domain, - doc="""Defining a substitute for conc_interface_MEA""", + doc="Volumetric liquid holdup [-]", ) - def sqrt_conc_interface_MEA_eqn(blk, t, x): + def holdup_liq_eqn(blk, t, x): if x == blk.liquid_phase.length_domain.last(): return Constraint.Skip else: - return ( - blk.sqrt_conc_interface_MEA[t, x] - == blk.conc_interface_MEA[t, x] ** 0.5 + log_holdup_parAlpha = log( + blk.holdup_parAlpha + * lunits("length") + * (lunits("acceleration") ** (2 / 3)) + ) + return blk.log_holdup_liq[t, x] == ( + blk.log_holdup_parA + + blk.holdup_parB + * ( + log_holdup_parAlpha + + blk.log_velocity_liq[t, x] + + (1 / 3) + * (blk.log_visc_d_liq[t, x] - blk.log_dens_mass_liq[t, x]) + ) ) + # ##################################################################### + + # Mass transfer coefficients of diffusing components in vapor phase [mol/m2.s.Pa] + # Mass transfer coefficients, Billet and Schultes (1999) correlation, + # where parameters are regressed by Chinen et al. (2018). + + self.Cv_ref = Var( + initialize=0.357, + units=pyunits.dimensionless, + doc="""Vapor packing specific constant in Billet and Schultes + volumetric mass transfer coefficient correlation""", + ) + + self.log_mass_transfer_coeff_vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + equilibrium_comp, + initialize=-10, + units=pyunits.dimensionless, + doc="Logarithm of the vapor mass transfer coeff", + ) + @self.Constraint( self.flowsheet().time, - self.liquid_phase.length_domain, - doc="""Enhancement factor - function of Hatta number""", + self.vapor_phase.length_domain, + equilibrium_comp, + doc="Defines log variable for the vapor mass transfer coeff", ) - def enhancement_factor_eqn1(blk, t, x): - if x == blk.liquid_phase.length_domain.last(): + def log_mass_transfer_coeff_vap_eqn(blk, t, x, j): + if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: - return blk.enhancement_factor[t, x] * ( - 1 - blk.conc_CO2_bulk[t, x] - ) == blk.Hatta[t, x] * (blk.sqrt_conc_interface_MEA[t, x]) * ( - 1 - blk.conc_CO2_equil_bulk[t, x] + return ( + exp(blk.log_mass_transfer_coeff_vap[t, x, j]) + * lunits("amount") + / lunits("pressure") + / lunits("length") ** 2 + / lunits("time") + == blk.mass_transfer_coeff_vap[t, x, j] ) + self.log_holdup_vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=-1, + units=pyunits.dimensionless, + doc="""Logarithm of vapor holdup""", + ) + @self.Constraint( self.flowsheet().time, - self.liquid_phase.length_domain, - doc="""Enhancement factor - function of instantaneous enhancement factor""", + self.vapor_phase.length_domain, + doc="""Defines log variable for vapor holdup""", ) - def enhancement_factor_eqn2(blk, t, x): - if x == blk.liquid_phase.length_domain.last(): + def log_holdup_vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: - return (blk.enhancement_factor[t, x] - 1) * ( - 1 - blk.conc_CO2_bulk[t, x] - ) == (blk.instant_E[t, x] - 1) * ( - 1 - - blk.sqrt_conc_interface_MEA[t, x] - * blk.sqrt_conc_interface_MEA[t, x] - ) + return exp(blk.log_holdup_vap[t, x]) == blk.holdup_vap[t, x] - @self.Objective() - def enhancement_factor_obj(blk): - time_set = self.flowsheet().time - x_set = blk.liquid_phase.length_domain - return sum( - sum( - ( - (blk.enhancement_factor[t, x] - 1) - * (1 - blk.conc_CO2_bulk[t, x]) - - (blk.instant_E[t, x] - 1) - * ( - 1 - - blk.sqrt_conc_interface_MEA[t, x] - * blk.sqrt_conc_interface_MEA[t, x] - ) - ) - ** 2 - for t in time_set - ) - for x in x_set - if x != blk.liquid_phase.length_domain.last() - ) + self.log_Cv_ref = Var( + bounds=(None, 100), + initialize=-1.030019497202498, + doc="Logarithm of the vapor mass transfer coeff param", + units=pyunits.dimensionless, + ) - self.enhancement_factor_obj.deactivate() - # Note: the objective function is only activated in the - # initialization routine in an intermediate step, to improve - # model convergence. It is not included in the unit model. + @self.Constraint( + doc="Defines log variable for the vapor mass transfer coeff param" + ) + def log_Cv_ref_eqn(blk): + return exp(blk.log_Cv_ref) == blk.Cv_ref - # Heat transfer coefficients, Chilton Colburn analogy - # Vapor-liquid heat transfer coefficient [J/m2.s.K] + self.log_parameter_var_eqn_map[self.log_Cv_ref] = self.log_Cv_ref_eqn - def rule_heat_transfer_coeff(blk, t, x): + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + equilibrium_comp, + doc="Vapor phase mass transfer coefficient", + ) + def mass_transfer_coeff_vap_eqn(blk, t, x, j): if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip else: - return blk.heat_transfer_coeff_base[ - t, x - ] == blk.mass_transfer_coeff_vap[ - t, x, "CO2" - ] * blk.vapor_phase.properties[ - t, x - ].pressure * blk.vapor_phase.properties[ - t, x - ].cp_mol_phase[ - "Vap" - ] * ( - blk.vapor_phase.properties[t, x].therm_cond_phase["Vap"] - / ( - blk.vapor_phase.properties[t, x].dens_mol_phase["Vap"] - * blk.vapor_phase.properties[t, x].cp_mol_phase["Vap"] - * blk.vapor_phase.properties[t, x].diffus_phase_comp[ - "Vap", "CO2" - ] + log_R_gas = log( + pyunits.convert(CONST.gas_constant, lunits("gas_constant")) + / lunits("gas_constant") + ) + log_vapor_temp = log( + pyunits.convert( + blk.vapor_phase.properties[t, x].temperature, + to_units=lunits("temperature"), + ) + / lunits("temperature") + ) + log_packing_specific_area = log( + blk.packing_specific_area + / (lunits("length") ** 2 / lunits("length") ** 3) + ) + log_hydraulic_diameter = log(blk.hydraulic_diameter / lunits("length")) + return blk.log_mass_transfer_coeff_vap[t, x, j] == ( + blk.log_Cv_ref + - log_R_gas + - log_vapor_temp + - 0.5 * blk.log_holdup_vap[t, x] + + 0.5 * (log_packing_specific_area - log_hydraulic_diameter) + + (2 / 3) * blk.log_diffus_vap_comp[t, x, j] + + (1 / 3) * (blk.log_visc_d_vap[t, x] - blk.log_dens_mass_vap[t, x]) + + (3 / 4) + * ( + blk.log_velocity_vap[t, x] + + blk.log_dens_mass_vap[t, x] + - log_packing_specific_area + - blk.log_visc_d_vap[t, x] ) - ) ** ( - 2 / 3 ) - self.heat_transfer_coeff_base_constraint = Constraint( - self.flowsheet().time, - self.vapor_phase.length_domain, - rule=rule_heat_transfer_coeff, - doc="""vap-liq heat transfer coefficient""", - ) + # ##################################################################### - # Flood point calculations - def rule_flood_velocity(blk, t, x): + # Mass transfer coefficients of diffusing components in liquid phase [m/s] + # Mass transfer coefficients, Billet and Schultes (1999) correlation, + # where parameters are regressed by Chinen et al. (2018). - if x == blk.vapor_phase.length_domain.first(): - return Expression.Skip + self.Cl_ref = Var( + initialize=0.5, + units=pyunits.dimensionless, + doc="""Liquid packing specific constant in Billet and Schultes + volumetric mass transfer coefficient correlation""", + ) - else: - x_liq = blk.vapor_phase.length_domain.prev(x) + self.log_mass_transfer_coeff_liq = Var( + self.flowsheet().time, + self.liquid_phase.length_domain, + solute_comp_list, + initialize=-9, + units=pyunits.dimensionless, + doc="""Logarithm of the liquid mass transfer coeff""", + ) - dens_liq = ( - blk.liquid_phase.properties[t, x_liq].mw - / blk.liquid_phase.properties[t, x_liq].vol_mol_phase["Liq"] + @self.Constraint( + self.flowsheet().time, + self.liquid_phase.length_domain, + solute_comp_list, + doc="""Defines log variable for the liquid mass transfer coeff""", + ) + def log_mass_transfer_coeff_liq_eqn(blk, t, x, j): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + exp(blk.log_mass_transfer_coeff_liq[t, x, j]) * lunits("velocity") + == blk.mass_transfer_coeff_liq[t, x, j] ) - dens_vap = ( - blk.vapor_phase.properties[t, x].mw - / blk.vapor_phase.properties[t, x].vol_mol_phase["Vap"] - ) + self.log_Cl_ref = Var( + bounds=(None, 100), + initialize=-0.30102999566, + units=pyunits.dimensionless, + doc="""Logarithm of the liquid mass transfer coeff param""", + ) - H = ( - blk.liquid_phase.properties[t, x_liq].flow_mass_phase["Liq"] - / blk.vapor_phase.properties[t, x].flow_mass_phase["Vap"] - ) * (dens_vap / dens_liq) ** 0.5 + @self.Constraint( + doc="""Defines log variable for the liquid mass transfer coeff param""" + ) + def log_Cl_ref_eqn(blk): + return exp(blk.log_Cl_ref) == blk.Cl_ref - mu_water = 0.001 * pyunits.Pa * pyunits.s + self.log_parameter_var_eqn_map[self.log_Cl_ref] = self.log_Cl_ref_eqn - return ( - ( - CONST.acceleration_gravity - * (blk.eps_ref) ** 3 - / blk.packing_specific_area - ) - * (dens_liq / dens_vap) + @self.Constraint( + self.flowsheet().time, + self.liquid_phase.length_domain, + solute_comp_list, + doc="Liquid phase mass transfer coefficient", + ) + def mass_transfer_coeff_liq_eqn(blk, t, x, j): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + log_hydraulic_diameter = log(blk.hydraulic_diameter / lunits("length")) + return blk.log_mass_transfer_coeff_liq[t, x, j] == ( + blk.log_Cl_ref + + (1 / 6) * log(12) + + 0.5 * ( - blk.liquid_phase.properties[t, x_liq].visc_d_phase["Liq"] - / mu_water + blk.log_velocity_liq[t, x] + + blk.log_diffus_liq_comp[t, x, j] + - blk.log_holdup_liq[t, x] + - log_hydraulic_diameter ) - ** (-0.2) - * exp(-4 * (H) ** 0.25) - ) ** 0.5 + ) + + # ##################################################################### + + # TODO We shouldn't be fixing variables in a construction method + # Fix mass transfer parameters + self.Cv_ref.fix() + self.Cl_ref.fix() - self.gas_velocity_fld = Expression( + # Fix interfacial area parameters + self.area_interfacial_parA.fix() + self.area_interfacial_parB.fix() + + # Fix liquid holdup parameters + self.holdup_parA.fix() + self.holdup_parB.fix() + + # --------------------------------------------------------------------- + # Vapor-liquid heat transfer coeff modified by Ackmann factor + self.heat_transfer_coeff_base = Var( self.flowsheet().time, self.vapor_phase.length_domain, - rule=rule_flood_velocity, - doc="Gas velocity at flooding point", + initialize=100, + units=lunits("power") / lunits("temperature") / lunits("length") ** 2, + doc="Uncorrected vapor-liquid heat transfer coefficient", ) - self.flood_fraction = Var( + self.log_heat_transfer_coeff_base = Var( self.flowsheet().time, self.vapor_phase.length_domain, - initialize=0.7, + initialize=5, units=pyunits.dimensionless, - doc="""Dimensionless flooding fraction""", ) - def rule_flood_fraction(blk, t, x): + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def log_heat_transfer_coeff_base_eqn(blk, t, x): if x == blk.vapor_phase.length_domain.first(): return Constraint.Skip - else: return ( - blk.flood_fraction[t, x] * blk.gas_velocity_fld[t, x] - == blk.velocity_vap[t, x] + exp(blk.log_heat_transfer_coeff_base[t, x]) + * lunits("heat_transfer_coefficient") + == blk.heat_transfer_coeff_base[t, x] ) - self.flood_fraction_constr = Constraint( + # Heat transfer coefficients, Chilton Colburn analogy + # Vapor-liquid heat transfer coefficient [J/m2.s.K] + @self.Constraint( self.flowsheet().time, self.vapor_phase.length_domain, - rule=rule_flood_fraction, - doc="Flooding fraction (expected to be less than 0.8)", + doc="""Define vap-liq heat transfer coefficient""", ) + def heat_transfer_coeff_base_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + return 3 * blk.log_heat_transfer_coeff_base[t, x] == 3 * ( + blk.log_mass_transfer_coeff_vap[t, x, "CO2"] + blk.log_pressure_vap[t, x] + ) + 2 * blk.log_therm_cond_vap[t, x] + blk.log_cp_mol_vap[t, x] - 2 * ( + blk.log_dens_mol_vap[t, x] + blk.log_diffus_vap_comp[t, x, "CO2"] + ) - # Scaling Routine - def calculate_scaling_factors_props(blk): - for x in blk.vapor_phase.length_domain: - iscale.set_scaling_factor(blk.vapor_phase.properties[0, x].pressure, 1e-5) - - iscale.set_scaling_factor(blk.vapor_phase.properties[0, x].flow_mol, 1e-4) - - for x in blk.liquid_phase.length_domain: - iscale.set_scaling_factor(blk.liquid_phase.properties[0, x].pressure, 1e-5) - - iscale.set_scaling_factor(blk.liquid_phase.properties[0, x].flow_mol, 1e-4) - - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].mole_frac_phase_comp_true[ - "Liq", "CO2" - ], - 1e4, - ) - - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].mole_frac_phase_comp_true[ - "Liq", "HCO3_-" - ], - 100, - ) - - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].mole_frac_phase_comp_true[ - "Liq", "MEACOO_-" - ], - 10, - ) - - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].mole_frac_phase_comp_true[ - "Liq", "MEA_+" - ], - 10, + if self.config.corrected_hx_coeff_eqn: + self.Ack_intermediate = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=1, + doc="Ackmann factor as intermediate variable", ) - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].mole_frac_phase_comp_true[ - "Liq", "MEA" - ], - 1, - ) + def rule_Ack_intermediate(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + return 1e-3 * blk.Ack_intermediate[ + t, x + ] * blk.heat_transfer_coeff_base[t, x] * blk.area_interfacial[ + t, x + ] * blk.area_column == -1e-3 * sum( + blk.vapor_phase.properties[t, x].cp_mol_phase_comp["Vap", j] + * blk.interphase_mass_transfer[t, x, j] + for j in equilibrium_comp + ) - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].mole_frac_phase_comp_true[ - "Liq", "H2O" - ], - 1, + self.Ack_intermediate_eqn = Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + rule=rule_Ack_intermediate, + doc="Constraint to calculate Ackmann factor", ) + # This formulation is bad. Because Ack_intermediate can be and often is near zero, the denominator is + # often near singular. Fortunately, in those cases, the numerator should also be near zero in such a + # way that the singularities cancel out, but that's only at a solved state, not intermediate iterations. + # Bringing (1 - exp(blk.Ack_intermediate[t, x])) to the other side just moves the singularity from the + # derivative calculation to the Jacobian inversion, which might allow IPOPT's regularization to save you + # but also causes a catastrophic loss of precision in the equation. + + # In the event we need to re-enable this equation, one can write h_corrected = h_base * Ack/(exp(Ack) - 1). + # When calculating x/(exp(x) - 1) directly, the numerator and denominator will both go to zero at the same + # rate at the same time, so should behave better more often. We can also take the approximation + # x/(exp(x) - 1) ~= 1 - x/2 + x**2/12 - x**4/720 + O(x**6) to avoid even worrying about division by zero. + + # Warning: blk.heat_transfer_coeff_base has units W/(m**2 K) as expected, but blk.heat_transfer_coeff has + # units W/(m K) because it isn't actually a heat transfer coefficient, but a heat transfer coefficient + # multiplied by heat transfer area per unit column length. + + # Please double check calculations before using them directly. + # Doug A. + def rule_heat_transfer_coeff(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + return blk.heat_transfer_coeff[t, x] == sum( + blk.vapor_phase.properties[t, x].cp_mol_phase_comp["Vap", j] + * blk.interphase_mass_transfer[t, x, j] + for j in equilibrium_comp + ) / (1 - exp(blk.Ack_intermediate[t, x])) - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].flow_mol_phase_comp_true[ - "Liq", "CO2" - ], - 100, - ) + else: - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].flow_mol_phase_comp_true[ - "Liq", "H2O" - ], - 1e-3, - ) + def rule_heat_transfer_coeff(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + return ( + blk.heat_transfer_coeff[t, x] + == blk.heat_transfer_coeff_base[t, x] + * blk.area_interfacial[t, x] + * blk.area_column + ) - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].flow_mol_phase_comp_true[ - "Liq", "MEA" - ], - 0.1, - ) + self.heat_transfer_coeff_eqn = Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + rule=rule_heat_transfer_coeff, + doc="Vap-Liq heat transfer correction by Ackmann factor", + ) - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].log_k_eq["carbamate"], 1 - ) - iscale.set_scaling_factor( - blk.liquid_phase.properties[0, x].log_k_eq["bicarbonate"], 1 + if self.config.surrogate_enhancement_factor_model is None: + self.enhancement_factor_vars, self.enhancement_factor_constraints = make_enhancement_factor_model( + self, + lunits, + kinetics="Luo" ) - - for v in blk.vapor_phase.properties.values(): - iscale.constraint_scaling_transform( - v.total_flow_balance, - iscale.get_scaling_factor(v.flow_mol, default=1, warning=True), + else: + self.CO2_loading = Var( + self.flowsheet().time, + self.liquid_phase.length_domain, + units=pyunits.dimensionless, + initialize=0.3 ) - - for v in blk.liquid_phase.properties.values(): - for (p, j) in v.appr_to_true_species.keys(): - iscale.constraint_scaling_transform( - v.appr_to_true_species[p, j], - iscale.get_scaling_factor( - v.flow_mol_phase_comp_true[p, j], default=1, warning=True - ), + @self.Constraint(self.flowsheet().time, self.liquid_phase.length_domain) + def CO2_loading_eqn(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip + return ( + self.CO2_loading[t, x] + * self.liquid_phase.properties[t, x].mole_frac_comp["MEA"] + == self.liquid_phase.properties[t, x].mole_frac_comp["CO2"] ) - for (p, j) in v.true_mole_frac_constraint.keys(): - iscale.constraint_scaling_transform( - v.true_mole_frac_constraint[p, j], - iscale.get_scaling_factor( - v.flow_mol_phase_comp_true[p, j], default=1, warning=True - ), + + self.H2O_loading = Var( + self.flowsheet().time, + self.liquid_phase.length_domain, + units=pyunits.dimensionless, + initialize=7.5 + ) + @self.Constraint(self.flowsheet().time, self.liquid_phase.length_domain) + def H2O_loading_eqn(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip + return ( + self.H2O_loading[t, x] + * self.liquid_phase.properties[t, x].mole_frac_comp["MEA"] + == self.liquid_phase.properties[t, x].mole_frac_comp["H2O"] ) - - iscale.constraint_scaling_transform(v.component_flow_balances["CO2"], 100) - iscale.constraint_scaling_transform(v.component_flow_balances["H2O"], 1) - iscale.constraint_scaling_transform(v.component_flow_balances["MEA"], 1) - - for v in blk.liquid_phase.properties.values(): - iscale.constraint_scaling_transform( - v.total_flow_balance, - iscale.get_scaling_factor(v.flow_mol, default=1, warning=True), + self.log_pCO2 = Var( + self.flowsheet().time, + self.liquid_phase.length_domain, # Going to phase shift this in the constraint + units=pyunits.dimensionless, + initialize=6, + doc="Logarithm of vapor phase vapor pressure of CO2" ) - - for v in blk.liquid_phase.properties.values(): - for rxn in v.log_k_eq_constraint.keys(): - iscale.constraint_scaling_transform( - v.log_k_eq_constraint[rxn], - iscale.get_scaling_factor(v.log_k_eq[rxn], default=1, warning=True), + @self.Constraint(self.flowsheet().time, self.liquid_phase.length_domain) + def log_pCO2_eqn(blk, t, x): + if x == blk.liquid_phase.length_domain.last(): + return Constraint.Skip + zf = blk.liquid_phase.length_domain.next(x) + pressure = pyunits.convert( + blk.vapor_phase.properties[t, zf].pressure, + to_units=lunits("pressure"), ) - - iscale.constraint_scaling_transform( - v.inherent_equilibrium_constraint[rxn], - iscale.get_scaling_factor(v.log_k_eq[rxn], default=1, warning=True), + return ( + exp(self.log_pCO2[t, x]) * lunits("pressure") + == pressure + * self.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] ) + self._temperature_liquid = Reference( + self.liquid_phase.properties[:,:].temperature, + ctype=Var + ) + self._mass_transfer_coeff_liq_CO2 = Reference( + self.mass_transfer_coeff_liq[:, :, "CO2"], + ctype=Var + ) + self._mass_transfer_coeff_vap_CO2 = Reference( + self.mass_transfer_coeff_vap[:, :, "CO2"], + ctype=Var + ) + + self.enhancement_factor_surrogate = SurrogateBlock( + self.flowsheet().time, + self.liquid_phase.length_domain, + concrete=True + ) + for t in self.flowsheet().time: + for x in self.liquid_phase.length_domain: + self.enhancement_factor_surrogate[t, x].build_model( + self.config.surrogate_enhancement_factor_model, + input_vars=[ + self.CO2_loading[t, x], + self.H2O_loading[t, x], + self.log_pCO2[t, x], + self._temperature_liquid[t, x], + self._mass_transfer_coeff_liq_CO2[t, x], + self._mass_transfer_coeff_vap_CO2[t, x] + ], + output_vars=[self.log_enhancement_factor[t, x]], + use_surrogate_bounds=False, + ) + if x == self.liquid_phase.length_domain.last(): + self.enhancement_factor_surrogate[t, x].deactivate() + self.enhancement_factor_vars = ComponentSet([ + self.CO2_loading, + self.H2O_loading, + self.log_pCO2, + self.log_enhancement_factor, + ]) + self.enhancement_factor_constraints = ComponentSet([ + self.CO2_loading_eqn, + self.H2O_loading_eqn, + self.log_pCO2_eqn, + self.enhancement_factor_surrogate, + ]) - def calculate_scaling_factors_control_vol(blk): - - # Scale control volume level variables - for x in blk.vapor_phase.length_domain: - iscale.set_scaling_factor(blk.vapor_phase.heat[0, x], 1e-6) + # Flood point calculations - iscale.set_scaling_factor(blk.vapor_phase.enthalpy_transfer[0, x], 1e-6) + self.log_flow_mass_Liq_Vap = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=1, + ) - iscale.set_scaling_factor(blk.vapor_phase._enthalpy_flow[0, x, "Vap"], 1e-6) + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Constraint for log variable", + ) + def log_flow_mass_Liq_Vap_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + x_liq = blk.liquid_phase.length_domain.prev(x) + return ( + exp(blk.log_flow_mass_Liq_Vap[t, x]) + * blk.vapor_phase.properties[t, x].flow_mass_phase["Vap"] + == blk.liquid_phase.properties[t, x_liq].flow_mass_phase["Liq"] + ) - iscale.set_scaling_factor( - blk.vapor_phase.enthalpy_flow_dx[0, x, "Vap"], 1e-6 - ) + self.log_flood_H = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=-1, + units=pyunits.dimensionless, + ) - for x in blk.liquid_phase.length_domain: - iscale.set_scaling_factor( - blk.liquid_phase._enthalpy_flow[0, x, "Liq"], 1e-8 - ) + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def flood_H_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + else: + x_liq = blk.liquid_phase.length_domain.prev(x) + return blk.log_flood_H[t, x] == blk.log_flow_mass_Liq_Vap[ + t, x + ] + 0.5 * ( + blk.log_dens_mass_vap[t, x] - blk.log_dens_mass_liq[t, x_liq] + ) - iscale.set_scaling_factor( - blk.liquid_phase.enthalpy_flow_dx[0, x, "Liq"], 1e-6 - ) + self.fourth_root_flood_H = Var( + self.flowsheet().time, self.vapor_phase.length_domain, initialize=1 + ) - iscale.set_scaling_factor(blk.liquid_phase.enthalpy_transfer[0, x], 1e-6) + @self.Constraint(self.flowsheet().time, self.vapor_phase.length_domain) + def fourth_root_flood_H_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip + return blk.fourth_root_flood_H[t, x] == exp(0.25 * blk.log_flood_H[t, x]) - iscale.set_scaling_factor(blk.liquid_phase.heat[0, x], 1e-6) + self.gas_velocity_fld = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + initialize=1, + units=lunits("velocity"), + doc="Flooding velocity", + ) - for (t, x, p, j), v in blk.vapor_phase._flow_terms.items(): - iscale.set_scaling_factor(blk.vapor_phase._flow_terms[t, x, p, j], 0.01) + self.log_gas_velocity_fld = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + bounds=(None, 100), + initialize=0.01, + ) - for (t, x, p, j), v in blk.vapor_phase.material_flow_dx.items(): - if x != 0: - iscale.set_scaling_factor( - blk.vapor_phase.material_flow_dx[t, x, p, j], 0.01 - ) + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Constraint for log variable", + ) + def log_gas_velocity_fld_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip else: - pass - - for (t, x, p, j), v in blk.liquid_phase._flow_terms.items(): - iscale.set_scaling_factor(blk.liquid_phase._flow_terms[t, x, p, j], 0.01) + return ( + exp(blk.log_gas_velocity_fld[t, x]) * lunits("velocity") + == blk.gas_velocity_fld[t, x] + ) - for (t, x, p, j), v in blk.liquid_phase.material_flow_dx.items(): - if x != 1: - if j != "MEA": - iscale.set_scaling_factor( - blk.liquid_phase.material_flow_dx[t, x, p, j], 0.01 - ) + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + ) + def gas_velocity_fld_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Expression.Skip else: - pass + x_liq = blk.liquid_phase.length_domain.prev(x) + log_packing_specific_area = log( + blk.packing_specific_area + / (lunits("length") ** 2 / lunits("length") ** 3) + ) + log_g = log( + pyunits.convert( + CONST.acceleration_gravity, to_units=lunits("acceleration") + ) + / lunits("acceleration") + ) + # Value of 0.001 Pa*s for viscosity of water was here when I found it: Doug + log_visc_d_liq_H2O = log( + pyunits.convert( + 0.001 * pyunits.Pa * pyunits.s, + to_units=lunits("dynamic_viscosity"), + ) + / lunits("dynamic_viscosity") + ) + return blk.log_gas_velocity_fld[t, x] == 0.5 * ( + log_g + + 3 * log(blk.eps_ref) + - log_packing_specific_area + + blk.log_dens_mass_liq[t, x_liq] + - blk.log_dens_mass_vap[t, x] + - 0.2 * (blk.log_visc_d_liq[t, x_liq] - log_visc_d_liq_H2O) + - 4 * blk.fourth_root_flood_H[t, x] + ) - for ( - t, - x, - p, - j, - ), v in blk.vapor_phase.material_flow_linking_constraints.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.vapor_phase._flow_terms[t, x, p, j], default=1, warning=True - ), - ) + self.flood_fraction = Var( + self.flowsheet().time, + self.vapor_phase.length_domain, + initialize=0.7, + units=pyunits.dimensionless, + doc="Dimensionless flooding fraction", + ) - for (t, x, j), v in blk.vapor_phase.material_balances.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.vapor_phase.material_flow_dx[t, x, "Vap", j], - default=1, - warning=True, - ), - ) + @self.Constraint( + self.flowsheet().time, + self.vapor_phase.length_domain, + doc="Flooding fraction (expected to be less than 0.8)", + ) + def flood_fraction_eqn(blk, t, x): + if x == blk.vapor_phase.length_domain.first(): + return Constraint.Skip - for (t, x, p), v in blk.vapor_phase.enthalpy_flow_linking_constraint.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.vapor_phase._enthalpy_flow[t, x, p], default=1, warning=True - ), - ) + else: + return ( + blk.flood_fraction[t, x] * blk.gas_velocity_fld[t, x] + == blk.velocity_vap[t, x] + ) - for (t, x), v in blk.vapor_phase.enthalpy_balances.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.vapor_phase.enthalpy_flow_dx[t, x, "Vap"], - default=1, - warning=True, - ), - ) + # Scaling Routine + def calculate_scaling_factors(self): + super().calculate_scaling_factors() - for (t, x), v in blk.vapor_phase.pressure_balance.items(): - iscale.constraint_scaling_transform(v, 1) + vunits = ( + self.config.vapor_phase.property_package.get_metadata().get_derived_units + ) + lunits = ( + self.config.liquid_phase.property_package.get_metadata().get_derived_units + ) - for (t, x, p, j), v in blk.vapor_phase.material_flow_dx_disc_eq.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.vapor_phase.material_flow_dx[t, x, p, j], - default=1, - warning=True, - ), - ) + def gsf(var): + return iscale.get_scaling_factor(var, default=1, warning=True) - for (t, x, p), v in blk.vapor_phase.enthalpy_flow_dx_disc_eq.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.vapor_phase.enthalpy_flow_dx[t, x, p], default=1, warning=True - ), - ) + def ssf(var, s): + iscale.set_scaling_factor(var, s, overwrite=False) - for (t, x), v in blk.vapor_phase.pressure_dx_disc_eq.items(): - iscale.constraint_scaling_transform(v, 1) + def cst(con, s): + iscale.constraint_scaling_transform(con, s, overwrite=False) - for ( - t, - x, - p, - j, - ), v in blk.liquid_phase.material_flow_linking_constraints.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase._flow_terms[t, x, p, j], default=1, warning=True - ), - ) + for t in self.flowsheet().time: + for x_liq in self.liquid_phase.length_domain: + if x_liq == self.liquid_phase.length_domain.last(): + continue + x_vap = self.vapor_phase.length_domain.next(x_liq) + sf_flow_mol = gsf(self.liquid_phase.properties[t, x_liq].flow_mol) - for (t, x, j), v in blk.liquid_phase.material_balances.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase.material_flow_dx[t, x, "Liq", j], - default=1, - warning=True, - ), - ) + cst(self.velocity_liq_eqn[t, x_liq], sf_flow_mol) + # Decent initial guess for liquid velocity + sf_v = iscale.get_scaling_factor( + self.velocity_liq[t, x_liq], default=0.016, warning=False + ) + ssf(self.log_velocity_liq[t, x_liq], 1) + cst(self.log_velocity_liq_eqn[t, x_liq], sf_v) - for (t, x, p), v in blk.liquid_phase.enthalpy_flow_linking_constraint.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase._enthalpy_flow[t, x, p], default=1, warning=True - ), - ) + sf_A_interfacial = iscale.get_scaling_factor( + self.area_interfacial[t, x_vap], + default=1 / value(self.packing_specific_area), + warning=False, + ) + # In the (common) case where it didn't have anything set, set the default value + ssf(self.area_interfacial[t, x_vap], sf_A_interfacial) + cst(self.log_area_interfacial_eqn[t, x_vap], sf_A_interfacial) + ssf(self.log_area_interfacial[t, x_vap], 1) + + sf_rho_l = iscale.get_scaling_factor( + self.liquid_phase.properties[t, x_liq].dens_mass_phase["Liq"], + default=1e-3, + warning=False, + ) + cst(self.log_dens_mass_liq_eqn[t, x_liq], sf_rho_l) - for (t, x), v in blk.liquid_phase.enthalpy_balances.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase.enthalpy_flow_dx[t, x, "Liq"], + sf_rho_v = iscale.get_scaling_factor( + self.vapor_phase.properties[t, x_vap].dens_mass_phase["Vap"], default=1, - warning=True, - ), - ) + warning=False, + ) + cst(self.log_dens_mass_vap_eqn[t, x_vap], sf_rho_v) - for (t, x, p, j), v in blk.liquid_phase.material_flow_dx_disc_eq.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase.material_flow_dx[t, x, p, j], - default=1, - warning=True, - ), - ) + sf = gsf(self.vapor_phase.properties[t, x_vap].pressure) + cst(self.log_pressure_vap_eqn[t, x_vap], sf) - for (t, x, p), v in blk.liquid_phase.enthalpy_flow_dx_disc_eq.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase.enthalpy_flow_dx[t, x, p], default=1, warning=True - ), - ) + sf = iscale.get_scaling_factor( + self.vapor_phase.properties[t, x_vap].dens_mol_phase["Vap"], + default=1 / 50, + warning=False, + ) + cst(self.log_dens_mol_vap_eqn[t, x_vap], sf) - def calculate_scaling_factors_unit_model(blk): + sf = iscale.get_scaling_factor( + self.vapor_phase.properties[t, x_vap].therm_cond_phase["Vap"], + default=100, + warning=False, + ) + cst(self.log_therm_cond_vap_eqn[t, x_vap], sf) - # Scale unit model level variables - for (t, x, j), v in blk.pressure_equil.items(): - if x != 0: - iscale.set_scaling_factor( - v, - 1 - / value(blk.liquid_phase.properties[t, x].fug_phase_comp["Liq", j]), + sf = iscale.get_scaling_factor( + self.vapor_phase.properties[t, x_vap].cp_mol_phase["Vap"], + default=1 / 50, + warning=False, ) - else: - iscale.set_scaling_factor(v, 1) + cst(self.log_cp_mol_vap_eqn[t, x_vap], sf) - for (t, x, j), v in blk.interphase_mass_transfer.items(): - if x != 0: - iscale.set_scaling_factor( - v, 1 / value(blk.interphase_mass_transfer[t, x, j]) + sf = iscale.get_scaling_factor( + self.heat_transfer_coeff_base[t, x_vap], + default=1 / 200, + warning=False, ) - else: - iscale.set_scaling_factor(v, 1) + ssf(self.heat_transfer_coeff_base[t, x_vap], sf) + cst(self.log_heat_transfer_coeff_base_eqn[t, x_vap], sf) + + for j in self.log_diffus_liq_comp_list: + sf_diffus_liq_comp = iscale.get_scaling_factor( + self.liquid_phase.properties[t, x_liq].diffus_phase_comp[ + "Liq", j + ], + default=2e8, + warning=False, + ) + cst(self.log_diffus_liq_comp_eqn[t, x_liq, j], sf_diffus_liq_comp) - # --------------------------------------------------------------------- - # Scale constraints + for j in self.solute_comp_list: + sf = iscale.get_scaling_factor( + self.mass_transfer_coeff_liq[t, x_liq, j], + default=1e4, + warning=False, + ) + ssf(self.mass_transfer_coeff_liq[t, x_liq, j], sf) + cst(self.log_mass_transfer_coeff_liq_eqn[t, x_liq, j], sf) - for (t, x), v in blk.mechanical_equilibrium.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.liquid_phase.properties[t, x].pressure, default=1, warning=False - ), - ) + sf = gsf(self.liquid_phase.properties[t, x_liq].visc_d_phase["Liq"]) + cst(self.log_visc_d_liq_eqn[t, x_liq], sf) - for (t, x, j), v in blk.pressure_at_interface.items(): - iscale.constraint_scaling_transform( - v, - iscale.get_scaling_factor( - blk.pressure_equil[t, x, j], default=1, warning=False - ), - ) + sf = gsf(self.vapor_phase.properties[t, x_vap].visc_d_phase["Vap"]) + cst(self.log_visc_d_vap_eqn[t, x_vap], sf) - for (t, x), v in blk.heat_transfer_eqn1.items(): + sf = gsf(self.vapor_phase.properties[t, x_vap].flow_mol) + cst(self.velocity_vap_eqn[t, x_vap], sf) + + sf = gsf(self.liquid_phase.properties[t, x_liq].flow_mass_phase["Liq"]) + cst(self.log_flow_mass_Liq_Vap_eqn[t, x_vap], sf) + + for j in self.equilibirum_comp: + sf = iscale.get_scaling_factor( + self.mass_transfer_coeff_vap[t, x_vap, j], + default=25000, + warning=False, + ) + ssf(self.mass_transfer_coeff_vap[t, x_vap, j], sf) + cst(self.log_mass_transfer_coeff_vap_eqn[t, x_vap, j], sf) + + sf = iscale.get_scaling_factor( + self.vapor_phase.properties[t, x_vap].diffus_phase_comp[ + "Vap", j + ], + default=5e4, + warning=False, + ) + cst(self.log_diffus_vap_comp_eqn[t, x_vap, j], sf) + + if self.config.surrogate_enhancement_factor_model is not None: + sf_x_CO2 = gsf(self.liquid_phase.properties[t, x_liq].mole_frac_comp["CO2"]) + sf_x_MEA = gsf(self.liquid_phase.properties[t, x_liq].mole_frac_comp["MEA"]) + sf_x_H2O = gsf(self.liquid_phase.properties[t, x_liq].mole_frac_comp["H2O"]) + + ssf(self.CO2_loading[t, x_liq], sf_x_CO2/sf_x_MEA) + cst(self.CO2_loading_eqn[t, x_liq], sf_x_CO2) + + ssf(self.H2O_loading[t, x_liq], sf_x_H2O/sf_x_MEA) + cst(self.H2O_loading_eqn[t, x_liq], sf_x_H2O) + + sf_units = pyunits.convert_value( + 1, + from_units=1 / lunits("pressure"), + to_units=1 / vunits("pressure") + ) + + sf_P_vap = sf_units * gsf(self.vapor_phase.properties[t, x_vap].pressure) + sf_CO2_vap = gsf(self.vapor_phase.properties[t, x_vap].mole_frac_comp["CO2"]) + cst(self.log_pCO2_eqn[t, x_liq], sf_P_vap*sf_CO2_vap) + + # TODO bring this into new form later + for (t, x), con in self.heat_transfer_coeff_eqn.items(): + # TODO Figure out how to calculate this + ssf(self.heat_transfer_coeff[t, x], 3e-6) + # Formulated to have the same scale as the heat transfer coefficient + cst(con, gsf(self.heat_transfer_coeff[t, x])) + + for (t, x), v in self.heat_transfer_eqn1.items(): iscale.constraint_scaling_transform( v, iscale.get_scaling_factor( - blk.vapor_phase.heat[t, x], default=1, warning=True + self.vapor_phase.heat[t, x], default=1, warning=True ), ) - for (t, x), v in blk.heat_transfer_eqn2.items(): + for (t, x), v in self.heat_transfer_eqn2.items(): iscale.constraint_scaling_transform( v, iscale.get_scaling_factor( - blk.liquid_phase.heat[t, x], default=1, warning=True + self.liquid_phase.heat[t, x], default=1, warning=True ), ) - for (t, x), v in blk.enthalpy_transfer_eqn1.items(): + for (t, x), v in self.enthalpy_transfer_eqn1.items(): iscale.constraint_scaling_transform( v, iscale.get_scaling_factor( - blk.vapor_phase.enthalpy_transfer[t, x], default=1, warning=True + self.vapor_phase.enthalpy_transfer[t, x], default=1, warning=True ), ) - for (t, x), v in blk.enthalpy_transfer_eqn2.items(): + for (t, x), v in self.enthalpy_transfer_eqn2.items(): iscale.constraint_scaling_transform( v, iscale.get_scaling_factor( - blk.liquid_phase.enthalpy_transfer[t, x], default=1, warning=True + self.liquid_phase.enthalpy_transfer[t, x], default=1, warning=True ), ) - def set_init_values_correlation_vars(blk, nfe): + + def set_init_values_correlation_vars(blk, nfe, mode): """ This method fixes the initial values of mass and heat transfer coefficients interfacial area, enhancement factor, and liquid holdup, based on @@ -1418,39 +1785,76 @@ def set_init_values_correlation_vars(blk, nfe): ] # Define enhancement factor known values - enhancement_factor_values = [ - 10.59208019, - 10.83945557, - 11.08225811, - 11.3273241, - 11.5773751, - 11.83479677, - 12.10215397, - 12.38234194, - 12.6787182, - 12.99526603, - 13.33681033, - 13.70931137, - 14.12027367, - 14.57932831, - 15.09908116, - 15.69637908, - 16.39424972, - 17.2249598, - 18.23499171, - 19.49343044, - 21.10665852, - 23.24523241, - 26.19540321, - 30.46301154, - 36.99441863, - 47.67140404, - 66.4626304, - 101.9815504, - 169.7954706, - 258.2290563, - 1, - ] + if mode == "absorber": + enhancement_factor_values = [ + 10.59208019, + 10.83945557, + 11.08225811, + 11.3273241, + 11.5773751, + 11.83479677, + 12.10215397, + 12.38234194, + 12.6787182, + 12.99526603, + 13.33681033, + 13.70931137, + 14.12027367, + 14.57932831, + 15.09908116, + 15.69637908, + 16.39424972, + 17.2249598, + 18.23499171, + 19.49343044, + 21.10665852, + 23.24523241, + 26.19540321, + 30.46301154, + 36.99441863, + 47.67140404, + 66.4626304, + 101.9815504, + 169.7954706, + 258.2290563, + 1, + ] + elif mode == "stripper": + enhancement_factor_values = [ + 23.68, + 22.39, + 21.25, + 20.24, + 19.78, + 18.93, + 18.15, + 17.45, + 17.12, + 16.80, + 16.21, + 15.93, + 15.41, + 15.16, + 14.92, + 14.47, + 14.25, + 14.05, + 13.85, + 13.66, + 13.48, + 13.29, + 13.12, + 12.8, + 12.64, + 12.5, + 12.36, + 12.25, + 12.21, + 12.38, + 1, + ] + else: + raise RuntimeError nfe_base = len(interfacial_area_values) - 1 x_nfe_list_base = [i / nfe_base for i in range(nfe_base + 1)] @@ -1510,11 +1914,13 @@ def interpolate_init_values(numdiscretepts, xdata, ydata): if x == blk.liquid_phase.length_domain.last(): blk.mass_transfer_coeff_liq[0, x, "CO2"].fix(0.001) blk.holdup_liq[0, x].fix(0.001) - blk.enhancement_factor[0, x].fix(1) + blk.log_enhancement_factor[0, x].fix(0) else: blk.mass_transfer_coeff_liq[0, x, "CO2"].fix(k_l_co2_values[i]) blk.holdup_liq[0, x].fix(0.01) - blk.enhancement_factor[0, x].fix(enhancement_factor_values[i]) + blk.log_enhancement_factor[0, x].fix( + log(enhancement_factor_values[i]) + ) # ========================================================================= # Model initialization routine @@ -1526,6 +1932,8 @@ def initialize( outlvl=idaeslog.NOTSET, solver=None, optarg=None, + mode="absorber", + monolithic_enhancement_factor_system=True, ): """ MEA solvent Packed Column initialization. @@ -1546,80 +1954,97 @@ def initialize( solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") # Set solver options - if optarg is None: optarg = {} - if "bound_push" not in optarg: - optarg["bound_push"] = 1e-10 - if "max_iter" not in optarg: - optarg["max_iter"] = 100 opt = get_solver(solver, optarg) - blk.calculate_scaling_factors_props() - blk.calculate_scaling_factors_control_vol() + initial_dof = degrees_of_freedom(blk) + + lunits = ( + blk.config.liquid_phase.property_package.get_metadata().get_derived_units + ) unit_constraints = [ "pressure_at_interface", "interphase_mass_transfer_eqn", "liquid_mass_transfer_eqn", "vapor_mass_transfer_eqn", - "heat_transfer_coeff_corr", + "Ack_intermediate_eqn", + "heat_transfer_coeff_eqn", "heat_transfer_eqn1", "heat_transfer_eqn2", "enthalpy_transfer_eqn1", "enthalpy_transfer_eqn2", ] - velocity_constraints = ["velocity_vap_eqn", "velocity_liq_eqn"] - - interfacial_area_constraint = ["area_interfacial_constraint"] - - liquid_holdup_constraint = ["holdup_liq_constraint"] + interfacial_area_constraints = [ + "log_area_interfacial_parA_eqn", + "area_interfacial_eqn", + "log_area_interfacial_eqn", + ] - mass_transfer_coeff_vap_constraints = ["mass_transfer_coeff_vap_constraint"] + liquid_holdup_constraints = [ + "holdup_liq_eqn", + "log_holdup_liq_eqn", + ] - mass_transfer_coeff_liq_constraints = ["mass_transfer_coeff_liq_constraint"] + mass_transfer_coeff_vap_constraints = [ + "log_mass_transfer_coeff_vap_eqn", + "log_holdup_vap_eqn", + "log_Cv_ref_eqn", + "mass_transfer_coeff_vap_eqn", + ] - enhancement_factor_constraints = [ - "sqrt_conc_interface_MEA_eqn", - "enhancement_factor_eqn1", - "enhancement_factor_eqn2", + mass_transfer_coeff_liq_constraints = [ + "log_mass_transfer_coeff_liq_eqn", + "log_Cl_ref_eqn", + "mass_transfer_coeff_liq_eqn", ] - heat_transfer_coefficient_constraint = ["heat_transfer_coeff_base_constraint"] + heat_transfer_coefficient_constraint = [ + "heat_transfer_coeff_base_eqn", + "log_heat_transfer_coeff_base_eqn", + ] - flooding_constraint = ["flood_fraction_constr"] + flooding_constraints = [ + "log_flood_H_eqn", + "log_flow_mass_Liq_Vap_eqn", + "flood_H_eqn", + "fourth_root_flood_H_eqn", + "log_gas_velocity_fld_eqn", + "gas_velocity_fld_eqn", + "flood_fraction_eqn", + ] # --------------------------------------------------------------------- # Deactivate unit model level constraints for c in blk.component_objects(Constraint, descend_into=True): if c.local_name in unit_constraints: c.deactivate() - if c.local_name in velocity_constraints: - c.deactivate() - if c.local_name in interfacial_area_constraint: + if c.local_name in interfacial_area_constraints: c.deactivate() - if c.local_name in liquid_holdup_constraint: + if c.local_name in liquid_holdup_constraints: c.deactivate() if c.local_name in mass_transfer_coeff_vap_constraints: c.deactivate() if c.local_name in mass_transfer_coeff_liq_constraints: c.deactivate() - if c.local_name in enhancement_factor_constraints: - c.deactivate() if c.local_name in heat_transfer_coefficient_constraint: c.deactivate() - if c.local_name in flooding_constraint: + if c.local_name in flooding_constraints: c.deactivate() + for constraint in blk.enhancement_factor_constraints: + constraint.deactivate() + # Fix variables nfe = blk.config.finite_elements - blk.set_init_values_correlation_vars(nfe) + blk.set_init_values_correlation_vars(nfe, mode) # Interface pressure - blk.pressure_equil.fix() + blk.mass_transfer_driving_force.fix() # Molar flux blk.interphase_mass_transfer.fix(0.0) @@ -1627,7 +2052,9 @@ def initialize( blk.liquid_phase.mass_transfer_term.fix(0.0) # Heat transfer rate - blk.heat_transfer_coeff.fix(0.0) + if blk.config.corrected_hx_coeff_eqn: + blk.Ack_intermediate.fix() + blk.heat_transfer_coeff.fix() blk.vapor_phase.heat.fix(0.0) blk.liquid_phase.heat.fix(0.0) blk.vapor_phase.enthalpy_transfer.fix(0.0) @@ -1637,40 +2064,33 @@ def initialize( blk.velocity_vap.fix() blk.velocity_liq.fix() + blk.area_interfacial_parA.fix() + blk.area_interfacial.fix() + blk.log_area_interfacial.fix() + + blk.holdup_parA.fix() + blk.holdup_liq.fix() + blk.log_holdup_liq.fix() + # Enhancement factor - blk.conc_interface_MEA.fix() - blk.sqrt_conc_interface_MEA.fix() + for var in blk.enhancement_factor_vars: + var.fix() # Flood fraction + blk.log_flood_H.fix() + blk.log_flow_mass_Liq_Vap.fix() + blk.fourth_root_flood_H.fix() + blk.gas_velocity_fld.fix() + blk.log_gas_velocity_fld.fix() blk.flood_fraction.fix() - # --------------------------------------------------------------------- - # Provide state arguments for property package initialization - vap_comp = blk.config.vapor_phase.property_package.component_list - liq_apparent_comp = [ - c[1] for c in blk.liquid_phase.properties.phase_component_set - ] + # Mass transfer coefficients + blk.Cl_ref.fix() + blk.log_mass_transfer_coeff_liq.fix() - if vapor_phase_state_args is None: - vapor_phase_state_args = { - "flow_mol": blk.vapor_inlet.flow_mol[0].value, - "temperature": blk.vapor_inlet.temperature[0].value, - "pressure": blk.vapor_inlet.pressure[0].value, - "mole_frac_comp": { - j: blk.vapor_inlet.mole_frac_comp[0, j].value for j in vap_comp - }, - } - - if liquid_phase_state_args is None: - liquid_phase_state_args = { - "flow_mol": blk.liquid_inlet.flow_mol[0].value, - "temperature": blk.liquid_inlet.temperature[0].value, - "pressure": blk.liquid_inlet.pressure[0].value, - "mole_frac_comp": { - j: blk.liquid_inlet.mole_frac_comp[0, j].value - for j in liq_apparent_comp - }, - } + blk.log_mass_transfer_coeff_vap.fix() + blk.log_holdup_vap.fix() + blk.Cv_ref.fix() init_log.info("Step 1: Property Package initialization") @@ -1692,25 +2112,59 @@ def initialize( hold_state=True, ) + # Calculate log variables of (supposed to be fixed) parameters + for var, eqn in blk.log_parameter_var_eqn_map.items(): + for k in eqn: + calculate_variable_from_constraint(var[k], eqn[k]) + + # Calculate these velocities first because their log counterparts are in log_properties_var_eqn_map + blk.velocity_liq.unfix() + blk.velocity_liq_eqn.activate() + for k in blk.velocity_liq_eqn: + calculate_variable_from_constraint( + blk.velocity_liq[k], blk.velocity_liq_eqn[k] + ) + + blk.velocity_vap.unfix() + blk.velocity_vap_eqn.activate() + for k in blk.velocity_vap_eqn: + calculate_variable_from_constraint( + blk.velocity_vap[k], blk.velocity_vap_eqn[k] + ) + + # Calculate log variables that can be directly calculated from properties from the property packages + # (and log variables for superficial velocity, which itself can be calculated from the initialized properties) + for var, eqn in blk.log_property_var_eqn_map.items(): + var.unfix() + eqn.activate() + for k in eqn: + calculate_variable_from_constraint(var[k], eqn[k]) + + calculate_variable_from_constraint( + blk.area_column, blk.column_cross_section_area_eqn + ) + + # --------------------------------------------------------------------- init_log.info("Step 2: Steady-State isothermal mass balance") with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(blk, tee=slc.tee) + assert_optimal_termination(res) init_log.info_high("Step 2: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- init_log.info("Step 3: Interface equilibrium") - # Activate interface pressure constraint - blk.pressure_equil.unfix() + blk.mass_transfer_driving_force.unfix() blk.pressure_at_interface.activate() with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(blk, tee=slc.tee) + assert_optimal_termination(res) init_log.info_high("Step 3 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- - init_log.info("Step 4a: Isothermal absoption") + init_log.info("Step 4a: Isothermal absorption") init_log.info_high("Calculating mass flux") # Unfix mass transfer terms @@ -1721,36 +2175,46 @@ def initialize( with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(blk, tee=slc.tee) - - init_log.info("Step 4b: Isothermal chemical absoption") + assert_optimal_termination(res) + init_log.info("Step 4b: Isothermal chemical absorption") init_log.info_high("Adding mass transfer to material balances") - blk.vapor_phase.mass_transfer_term.unfix() blk.liquid_phase.mass_transfer_term.unfix() - blk.vapor_mass_transfer_eqn.activate() blk.liquid_mass_transfer_eqn.activate() + blk.vapor_phase.mass_transfer_term.unfix() + blk.vapor_mass_transfer_eqn.activate() + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(blk, tee=slc.tee) + assert_optimal_termination(res) init_log.info_high("Step 4 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- - init_log.info("Step 5: Adiabatic chemical absoption") + init_log.info("Step 5: Adiabatic chemical absorption") init_log.info_high("Isothermal to Adiabatic ") # Unfix heat transfer terms + if blk.config.corrected_hx_coeff_eqn: + blk.Ack_intermediate.unfix() + blk.Ack_intermediate_eqn.activate() blk.heat_transfer_coeff.unfix() - for c in ["heat_transfer_coeff_corr"]: - getattr(blk, c).activate() + blk.heat_transfer_coeff_eqn.activate() - for k in blk.heat_transfer_coeff_corr: + if blk.config.corrected_hx_coeff_eqn: + for k in blk.Ack_intermediate_eqn: + calculate_variable_from_constraint( + blk.Ack_intermediate[k], blk.Ack_intermediate_eqn[k] + ) + + for k in blk.heat_transfer_coeff_eqn: calculate_variable_from_constraint( - blk.heat_transfer_coeff[k], blk.heat_transfer_coeff_corr[k] + blk.heat_transfer_coeff[k], blk.heat_transfer_coeff_eqn[k] ) - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) - + # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + # res = opt.solve(blk, tee=slc.tee) + # assert_optimal_termination(res) blk.vapor_phase.heat.unfix() blk.liquid_phase.heat.unfix() blk.vapor_phase.enthalpy_transfer.unfix() @@ -1767,24 +2231,47 @@ def initialize( with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(blk, tee=slc.tee) - + assert_optimal_termination(res) init_log.info_high("Step 5 complete: {}.".format(idaeslog.condition(res))) - # --------------------------------------------------------------------- - - blk.calculate_scaling_factors_unit_model() - iscale.calculate_scaling_factors(blk) + init_log.info("Step 6: Flooding velocity constraints") - # --------------------------------------------------------------------- - init_log.info("Step 6: Velocity constraints") - blk.velocity_vap.unfix() - blk.velocity_liq.unfix() + blk.log_flood_H.unfix() + blk.log_flow_mass_Liq_Vap.unfix() + blk.fourth_root_flood_H.unfix() + blk.gas_velocity_fld.unfix() + blk.log_gas_velocity_fld.unfix() + blk.flood_fraction.unfix() - for c in velocity_constraints: - getattr(blk, c).activate() + blk.flood_H_eqn.activate() + blk.log_flow_mass_Liq_Vap_eqn.activate() + blk.flood_H_eqn.activate() + blk.fourth_root_flood_H_eqn.activate() + blk.log_gas_velocity_fld_eqn.activate() + blk.gas_velocity_fld_eqn.activate() + blk.flood_fraction_eqn.activate() + + for t in blk.flowsheet().time: + for x in blk.vapor_phase.length_domain: + if x == blk.vapor_phase.length_domain.first(): + pass + else: + calculate_variable_from_constraint( + blk.log_flood_H[t, x], blk.flood_H_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_flow_mass_Liq_Vap[t, x], + blk.log_flow_mass_Liq_Vap_eqn[t, x], + ) + calculate_variable_from_constraint( + blk.log_gas_velocity_fld[t, x], + blk.log_gas_velocity_fld_eqn[t, x], + ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 6 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- init_log.info("Step 7: Interfacial area constraint") @@ -1793,23 +2280,22 @@ def initialize( degrees_of_freedom(blk) ) ) - blk.area_interfacial.unfix() - for c in interfacial_area_constraint: - getattr(blk, c).activate() + # Confusing naming convention: log_var_eqn are always of the form exp(log_var) == var. + # log_area_interfacial is defined from the performance equation area_interfacial_eqn + # and then area_interfacial is back-calculated from the exponential relationship + for k in blk.log_area_interfacial_eqn: + calculate_variable_from_constraint( + blk.log_area_interfacial[k], blk.area_interfacial_eqn[k] + ) + calculate_variable_from_constraint( + blk.area_interfacial[k], blk.log_area_interfacial_eqn[k] + ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) - - # Scale variable - for (t, x), v in blk.area_interfacial.items(): - iscale.set_scaling_factor(v, 1 / value(blk.area_interfacial[t, x])) - - # Scale constraint - for (t, x), v in blk.area_interfacial_constraint.items(): - iscale.constraint_scaling_transform( - v, iscale.get_scaling_factor(blk.area_interfacial[t, x]) - ) + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 7 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- init_log.info("Step 8: Liquid holdup constraint") @@ -1818,31 +2304,20 @@ def initialize( degrees_of_freedom(blk) ) ) - blk.holdup_liq.unfix() - - for c in liquid_holdup_constraint: - getattr(blk, c).activate() - for x in blk.liquid_phase.length_domain: - if x == blk.liquid_phase.length_domain.last(): - pass - else: - calculate_variable_from_constraint( - blk.holdup_liq[0, x], blk.holdup_liq_constraint[0, x] - ) + # Same thing with holdup_liq + for k in blk.log_holdup_liq_eqn: + calculate_variable_from_constraint( + blk.log_holdup_liq[k], blk.holdup_liq_eqn[k] + ) + calculate_variable_from_constraint( + blk.holdup_liq[k], blk.log_holdup_liq_eqn[k] + ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) - - # Scale variable - for (t, x), v in blk.holdup_liq.items(): - iscale.set_scaling_factor(v, 1 / value(blk.holdup_liq[t, x])) - - # Scale constraint - for (t, x), v in blk.holdup_liq_constraint.items(): - iscale.constraint_scaling_transform( - v, iscale.get_scaling_factor(blk.holdup_liq[t, x]) - ) + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 8 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- init_log.info("Step 9: Vapor phase mass transfer coefficient constraint") @@ -1851,38 +2326,37 @@ def initialize( degrees_of_freedom(blk) ) ) + + blk.log_holdup_vap.unfix() + blk.log_mass_transfer_coeff_vap.unfix() blk.mass_transfer_coeff_vap.unfix() for c in mass_transfer_coeff_vap_constraints: getattr(blk, c).activate() - for x in blk.vapor_phase.length_domain: - if x == blk.vapor_phase.length_domain.first(): - pass - else: - calculate_variable_from_constraint( - blk.mass_transfer_coeff_vap[0, x, "CO2"], - blk.mass_transfer_coeff_vap_constraint[0, x, "CO2"], - ) - calculate_variable_from_constraint( - blk.mass_transfer_coeff_vap[0, x, "H2O"], - blk.mass_transfer_coeff_vap_constraint[0, x, "H2O"], - ) - - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) + for t in blk.flowsheet().time: + for x in blk.vapor_phase.length_domain: + if x == blk.liquid_phase.length_domain.first(): + pass + else: + calculate_variable_from_constraint( + blk.log_holdup_vap[t, x], blk.log_holdup_vap_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_mass_transfer_coeff_vap[t, x, "CO2"], + blk.mass_transfer_coeff_vap_eqn[t, x, "CO2"], + ) + calculate_variable_from_constraint( + blk.mass_transfer_coeff_vap[t, x, "H2O"], + blk.log_mass_transfer_coeff_vap_eqn[t, x, "H2O"], + ) - # Scale variable - for (t, x, j), v in blk.mass_transfer_coeff_vap.items(): - iscale.set_scaling_factor( - v, 1 / value(blk.mass_transfer_coeff_vap[t, x, j]) - ) + calculate_variable_from_constraint(blk.Cv_ref, blk.log_Cv_ref_eqn) - # Scale constraint - for (t, x, j), v in blk.mass_transfer_coeff_vap_constraint.items(): - iscale.constraint_scaling_transform( - v, iscale.get_scaling_factor(blk.mass_transfer_coeff_vap[t, x, j]) - ) + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 9 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- init_log.info("Step 10: Liquid phase mass transfer coefficient constraint") @@ -1893,141 +2367,144 @@ def initialize( ) blk.mass_transfer_coeff_liq.unfix() + # blk.log_diffus_liq_comp.unfix() + blk.log_mass_transfer_coeff_liq.unfix() + # blk.Cl_ref.unfix() + blk.mass_transfer_coeff_liq.unfix() + for c in mass_transfer_coeff_liq_constraints: getattr(blk, c).activate() - for x in blk.liquid_phase.length_domain: - if x == blk.liquid_phase.length_domain.last(): - pass - else: - calculate_variable_from_constraint( - blk.mass_transfer_coeff_liq[0, x, "CO2"], - blk.mass_transfer_coeff_liq_constraint[0, x, "CO2"], - ) + for t in blk.flowsheet().time: + for x in blk.liquid_phase.length_domain: + if x == blk.liquid_phase.length_domain.last(): + pass + else: + calculate_variable_from_constraint( + blk.log_mass_transfer_coeff_liq[t, x, "CO2"], + blk.mass_transfer_coeff_liq_eqn[t, x, "CO2"], + ) + calculate_variable_from_constraint( + blk.mass_transfer_coeff_liq[t, x, "CO2"], + blk.log_mass_transfer_coeff_liq_eqn[t, x, "CO2"], + ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) - - # Scale variable - for (t, x, j), v in blk.mass_transfer_coeff_liq.items(): - iscale.set_scaling_factor( - v, 1 / value(blk.mass_transfer_coeff_liq[t, x, "CO2"]) - ) - - # Scale constraint - for (t, x, j), v in blk.mass_transfer_coeff_liq_constraint.items(): - iscale.constraint_scaling_transform( - v, iscale.get_scaling_factor(blk.mass_transfer_coeff_liq[t, x, "CO2"]) - ) + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 10 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- - init_log.info("Step 11: Enhancement factor model constraints") init_log.info( - "Initializing enhancement factor (liquid phase) - degrees_of_freedom = {}".format( + "Step 11a Initializing enhancement factor model - degrees_of_freedom = {}".format( degrees_of_freedom(blk) ) ) - blk.conc_interface_MEA.unfix() - blk.sqrt_conc_interface_MEA.unfix() - blk.enhancement_factor.unfix() - - for c in enhancement_factor_constraints: - getattr(blk, c).activate() - - blk.enhancement_factor_eqn2.deactivate() - blk.enhancement_factor_obj.activate() - for (t, x), v in blk.enhancement_factor_eqn1.items(): - if x != 1: - iscale.constraint_scaling_transform(v, 100) - else: - pass - - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) - - blk.enhancement_factor_eqn2.activate() - - for (t, x), v in blk.enhancement_factor_eqn2.items(): - if x != 1: - iscale.constraint_scaling_transform(v, 100) - else: - pass + for constraint in blk.enhancement_factor_constraints: + constraint.activate() + for var in blk.enhancement_factor_vars: + var.unfix() + + if blk.config.surrogate_enhancement_factor_model is None: + initialize_enhancement_factor_model( + blk, + outlvl=outlvl, + optarg=optarg, + solver=solver, + ) + else: + for t in blk.flowsheet().time: + for x in blk.liquid_phase.length_domain: + if x == blk.liquid_phase.length_domain.last(): + blk.enhancement_factor_surrogate[t, x].deactivate() + continue + zf = blk.liquid_phase.length_domain.next(x) + blk.CO2_loading[t, x].set_value( + value( + blk.liquid_phase.properties[t, x].mole_frac_comp["CO2"] + / blk.liquid_phase.properties[t, x].mole_frac_comp["MEA"] + ) + ) + blk.H2O_loading[t, x].set_value( + value( + blk.liquid_phase.properties[t, x].mole_frac_comp["H2O"] + / blk.liquid_phase.properties[t, x].mole_frac_comp["MEA"] + ) + ) + pressure = pyunits.convert( + blk.vapor_phase.properties[t, zf].pressure, + to_units=lunits("pressure") + ) + blk.log_pCO2[t, x].set_value( + value( + log( + pressure + * blk.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + ) + ) + ) + + - blk.enhancement_factor_obj.deactivate() + init_log.info("Step 11b: Solve model with initialized enhancement factor") with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 11 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- init_log.info("Step 12: Heat transfer coefficient constraint") init_log.info( - "Initializing heat transfer coefficient - degrees_of_freedom = {}".format( + "Initializing heat transfer coefficent - degrees_of_freedom = {}".format( degrees_of_freedom(blk) ) ) blk.heat_transfer_coeff_base.unfix() - blk.heat_transfer_coeff_base_constraint.activate() + blk.log_heat_transfer_coeff_base.unfix() + blk.heat_transfer_coeff_base_eqn.activate() + blk.log_heat_transfer_coeff_base_eqn.activate() for x in blk.vapor_phase.length_domain: if x == blk.vapor_phase.length_domain.first(): pass else: + calculate_variable_from_constraint( + blk.log_heat_transfer_coeff_base[0, x], + blk.heat_transfer_coeff_base_eqn[0, x], + ) calculate_variable_from_constraint( blk.heat_transfer_coeff_base[0, x], - blk.heat_transfer_coeff_base_constraint[0, x], + blk.log_heat_transfer_coeff_base_eqn[0, x], ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) - - # Scale variable - for (t, x), v in blk.heat_transfer_coeff_base.items(): - iscale.set_scaling_factor(v, 1 / value(blk.heat_transfer_coeff_base[t, x])) - - # Scale constraint - for (t, x), v in blk.heat_transfer_coeff_base_constraint.items(): - iscale.constraint_scaling_transform( - v, iscale.get_scaling_factor(blk.heat_transfer_coeff_base[t, x]) - ) - - # --------------------------------------------------------------------- - init_log.info("Step 13: Flooding constraint") - init_log.info( - "Initializing flooding calculations - degrees_of_freedom = {}".format( - degrees_of_freedom(blk) - ) - ) - - blk.flood_fraction.unfix() - blk.flood_fraction_constr.activate() - - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(blk, tee=slc.tee) + res = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + init_log.info_high("Step 12 complete: {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- # Release state blk.vapor_phase.release_state(flags=vflag) blk.liquid_phase.release_state(flags=lflag) - init_log.info( - "Step 14: Initialization completed - degrees_of_freedom = {}".format( + init_log.info_high( + "Initialization completed - degrees_of_freedom = {}".format( degrees_of_freedom(blk) ) ) # Check DOF - if degrees_of_freedom(blk) != 0: + if degrees_of_freedom(blk) != initial_dof: raise InitializationError( - f"Degrees of freedom were not 0 at the end " + f"Degrees of freedom were not returned to their initial value of {initial_dof} at the end " f"of initialization. DoF = {degrees_of_freedom(blk)}" ) # Check solver status if not check_optimal_termination(res): raise InitializationError( - "Model failed to initialize successfully. Please check " - "the output logs for more information." + f"Model failed to initialize successfully. Please check " + f"the output logs for more information." ) - else: - init_log.info("Initialization successful") diff --git a/idaes/models_extra/column_models/enhancement_factor_model.py b/idaes/models_extra/column_models/enhancement_factor_model.py new file mode 100644 index 0000000000..8fca2b78c8 --- /dev/null +++ b/idaes/models_extra/column_models/enhancement_factor_model.py @@ -0,0 +1,624 @@ +# Import Python libraries +import logging +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd + +# Import Pyomo libraries +import pyomo.opt +from pyomo.environ import ( + Block, + ConcreteModel, + value, + Var, + Reals, + NonNegativeReals, + Param, + TransformationFactory, + Constraint, + Expression, + Objective, + SolverStatus, + TerminationCondition, + check_optimal_termination, + assert_optimal_termination, + exp, + log, + sqrt, + units as pyunits, + Set, + Reference +) +from pyomo.common.collections import ComponentSet, ComponentMap + +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.common.config import ConfigValue, Bool + +# Import IDAES Libraries +from idaes.core.util.constants import Constants as CONST +from idaes.models_extra.column_models.solvent_column import PackedColumnData + +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.initialization import _fix_vars, _restore_fixedness +from idaes.core import declare_process_block_class, FlowsheetBlock, StateBlock +from idaes.core.util.exceptions import InitializationError +from idaes.core.solvers.get_solver import get_solver +import idaes.logger as idaeslog + +from idaes.core.solvers import use_idaes_solver_configuration_defaults +import idaes.core.util.scaling as iscale +from pyomo.util.subsystems import ( + create_subsystem_block, +) +from idaes.core.solvers.petsc import ( + _sub_problem_scaling_suffix, +) + +def make_enhancement_factor_model(blk, lunits): + """ + Enhancement factor based liquid phase mass transfer model. + """ + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Second order rate constant [m3/(mol.s)]", + ) + def log_rate_constant(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + T = pyunits.convert( + b.liquid_phase.properties[t, x].temperature, to_units=pyunits.K + ) + C_MEA = pyunits.convert( + b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ], + to_units=pyunits.mol / pyunits.m ** 3, + ) + + # Putta, Svendsen, Knuutila 2017 Eqn. 42 + return log( + pyunits.convert( + 3.1732e9 + * exp(-4936.6 * pyunits.K / T) + * C_MEA + * 1e-6 + * ((pyunits.m) ** 6 / (pyunits.mol ** 2 * pyunits.s)), + to_units=1 / (lunits("time") * lunits("density_mole")), + ) + * lunits("time") + * lunits("density_mole") + ) + + blk.log_hatta_number = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=3, + doc="Hatta number", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Hatta number constraint", + ) + def hatta_number_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return b.log_hatta_number[t, x] == 0.5 * ( + b.log_rate_constant[t, x] + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + + b.log_diffus_liq_comp[t, x, "CO2"] + ) - b.log_mass_transfer_coeff_liq[t, x, "CO2"] + + blk.conc_CO2_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=0.95, + units=pyunits.dimensionless, + bounds=(0, None), + doc="""Dimensionless concentration of CO2, + Driving force term where + Absorption implies conc_CO2_bulk < 1 and + Desorption implies conc_CO2_bulk > 1 """, + ) + blk.log_conc_CO2_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=-0.25, + units=pyunits.dimensionless, + doc="Natural logarithm of conc_CO2_bulk", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Equation defining log_conc_CO2_bulk", + ) + def log_conc_CO2_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return exp(b.log_conc_CO2_bulk[t, x]) == b.conc_CO2_bulk[t, x] + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Dimensionless concentration of CO2""", + ) + def conc_CO2_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + zf = b.liquid_phase.length_domain.next(x) + Pressure = pyunits.convert( + b.vapor_phase.properties[t, zf].pressure, + to_units=lunits("pressure"), + ) + return b.log_conc_CO2_bulk[t, x] + log( + ( + b.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + * Pressure + / b.psi[t, zf] + + b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + ) / lunits("density_mole") + ) == b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + log( + b.liquid_phase.properties[t, x].henry["Liq", "CO2"] + / b.psi[t, zf] + + 1 + ) + # return ( + # b.conc_CO2_bulk[t, x] * ( + # b.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + # * Pressure + # / b.psi[t, zf] + # + b.liquid_phase.properties[t, x].conc_mol_phase_comp_true["Liq", "CO2"] + # ) + # == b.liquid_phase.properties[t, x].conc_mol_phase_comp_true["Liq", "CO2"] + # * ( + # b.liquid_phase.properties[t, x].henry["Liq", "CO2"] + # / b.psi[t, zf] + # + 1 + # ) + # ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Instantaneous Enhancement factor""", + ) + def log_instant_E_minus_one(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + zf = b.liquid_phase.length_domain.next(x) + return ( + b.log_diffus_liq_comp[t, x, "MEA"] + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + + b.log_conc_CO2_bulk[t, x] + - log(2) + - b.log_diffus_liq_comp[t, x, "CO2"] + - b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + ) + # ====================================================================== + # Enhancement factor model + # Reference: Jozsef Gaspar,Philip Loldrup Fosbol, (2015) + + blk.conc_interface_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + bounds=(0, None), + initialize=0.95, + units=pyunits.dimensionless, + doc="""Dimensionless concentration of MEA + at interface """, + ) + + blk.log_conc_interface_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + bounds=(None, 100), + initialize=-0.05, + units=pyunits.dimensionless, + doc="""Logarithm of conc_interface_MEA""", + ) + # ============================================================================= + # ------------------------ ORIGINAL ----------------------------------- + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Dimensionless concentration of MEACOO-", + ) + def conc_interface_MEACOO(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return 1 + ( + b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + ) * (1 - b.conc_interface_MEA[t, x]) / ( + 2 + * b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEACOO_-" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEACOO_-" + ] + ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Dimensionless concentration of MEA+", + ) + def conc_interface_MEAH(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return 1 + ( + b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + ) * (1 - b.conc_interface_MEA[t, x]) / ( + 2 + * b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA_+" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA_+" + ] + ) + + # ============================================================================= + blk.conc_CO2_equil_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=1, + bounds=(0, None), + doc="""Dimensionless concentration of CO2 + at equilibrium with the bulk """, + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Constraint for dimensionless concentration of CO2 + at equilibrium with the bulk """, + ) + def conc_CO2_equil_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + b.conc_CO2_equil_bulk[t, x] * b.conc_interface_MEA[t, x] ** 2 + == b.conc_CO2_bulk[t, x] + * b.conc_interface_MEAH[t, x] + * b.conc_interface_MEACOO[t, x] + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Defining log_conc_interface_MEA", + ) + def log_conc_interface_MEA_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + exp(b.log_conc_interface_MEA[t, x]) + == b.conc_interface_MEA[t, x] + ) + + blk.log_singular_CO2_CO2_ratio = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=-3, + bounds=(-12, 12), + doc="""Logarithm of the ratio of one minus dimensionless concentration CO2 at equilibrium + with the bulk to one minus dimensionless concentration of CO2 actually in the bulk. + Individually, numerator and denominator go to zero at equilibrium, but the ratio + remains (hopefully) well defined mathematically, if not numerically.""", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Equation defining log_singular_CO2_CO2_ratio", + ) + def log_singular_CO2_CO2_ratio_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + # Constraint written presently to lose meaning when dimensionless concentrations approach 1 + # If this form doesn't work, we can try one with division in it + return exp(b.log_singular_CO2_CO2_ratio[t, x]) * ( + 1 - b.conc_CO2_bulk[t, x] + ) == (1 - b.conc_CO2_equil_bulk[t, x]) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Enhancement factor - function of Hatta number", + ) + def enhancement_factor_eqn1(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + b.log_enhancement_factor[t, x] + == b.log_hatta_number[t, x] + + 0.5 * b.log_conc_interface_MEA[t, x] + + b.log_singular_CO2_CO2_ratio[t, x] + ) + # return (b.enhancement_factor[t, x] * ( + # 1 - b.conc_CO2_bulk[t, x] + # ) == b.Hatta[t, x] * b.sqrt_conc_interface_MEA[t, x] * ( + # 1 - b.conc_CO2_equil_bulk[t, x] + # ) + + # blk.log_singular_MEA_CO2_ratio = Var( + # blk.flowsheet().time, + # blk.liquid_phase.length_domain, + # initialize=-5, + # bounds=( + # -12, + # 3, + # ), # Should be straight up nonpositive, but give it a bit of slack + # doc="""Logarithm of the ratio of one minus dimensionless concentration MEA at interface + # to one minus dimensionless concentration of CO2 at equilibrium with the bulk. + # Individually, numerator and denominator go to zero at equilibrium, but the ratio + # remains (hopefully) well defined mathematically, if not numerically.""", + # ) + + # @blk.Constraint( + # blk.flowsheet().time, + # blk.liquid_phase.length_domain, + # doc="Equation defining log_singular_MEA_CO2_ratio", + # ) + # def log_singular_MEA_CO2_ratio_eqn(b, t, x): + # if x == b.liquid_phase.length_domain.last(): + # return Constraint.Skip + # else: + # # Constraint written presently to lose meaning when dimensionless concentrations approach 1 + # # If this form doesn't work, we can try one with division in it + # return exp(b.log_singular_MEA_CO2_ratio[t, x]) * ( + # 1 - b.conc_CO2_bulk[t, x] + # ) == (1 - b.conc_interface_MEA[t, x]) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Enhancement factor - function of instantaneous enhancement factor", + ) + def enhancement_factor_eqn2(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + # return b.log_enhancement_factor_minus_one[t, x] == ( + # b.log_instant_E_minus_one[t, x] + # + b.log_singular_MEA_CO2_ratio[t, x] + # ) + return (b.enhancement_factor[t, x] - 1) * ( + 1 - b.conc_CO2_bulk[t, x] + ) == exp(b.log_instant_E_minus_one[t, x]) * (1 - b.conc_interface_MEA[t, x]) + + enhancement_factor_vars = [ + blk.log_enhancement_factor, + blk.log_hatta_number, + blk.conc_CO2_bulk, + blk.log_conc_CO2_bulk, + blk.conc_interface_MEA, + blk.log_conc_interface_MEA, + blk.conc_CO2_equil_bulk, + blk.log_singular_CO2_CO2_ratio, + ] + + enhancement_factor_constraints = [ + blk.hatta_number_eqn, + blk.log_conc_CO2_bulk_eqn, + blk.conc_CO2_bulk_eqn, + blk.conc_CO2_equil_bulk_eqn, + blk.log_conc_interface_MEA_eqn, + blk.log_singular_CO2_CO2_ratio_eqn, + blk.enhancement_factor_eqn1, + blk.enhancement_factor_eqn2, + ] + + return enhancement_factor_vars, enhancement_factor_constraints + +def initialize_enhancement_factor_model( + blk, + state_args=None, + outlvl=idaeslog.NOTSET, + optarg=None, + solver=None, + ): + # Set up logger for initialization and solve + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") + # Set solver options + if optarg is None: + optarg = {} + + solver_obj = get_solver(solver, optarg) + + long_var_list = [] + long_eqn_list = [] + + for t in blk.flowsheet().time: + for x in blk.liquid_phase.length_domain: + if x == blk.liquid_phase.length_domain.last(): + continue + zf = blk.liquid_phase.length_domain.next(x) + calculate_variable_from_constraint( + blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_hatta_number[t, x], blk.hatta_number_eqn[t, x] + ) + + if value(blk.conc_CO2_bulk[t, x]) < 1: + blk.conc_interface_MEA[t, x].value = 0.95 + blk.log_conc_interface_MEA[t, x].value = log(0.95) + blk.log_singular_CO2_CO2_ratio[t, x].value = log(1.1) + else: + blk.conc_interface_MEA[t, x].value = 1.05 + blk.log_conc_interface_MEA[t, x].value = log(1.05) + blk.log_singular_CO2_CO2_ratio[t, x].value = log(0.9) + calculate_variable_from_constraint( + blk.conc_CO2_equil_bulk[t, x], + blk.conc_CO2_equil_bulk_eqn[t, x], + ) + + entangled_vars = [ + blk.conc_interface_MEA[t, x], + blk.log_conc_interface_MEA[t, x], + blk.log_enhancement_factor[t, x], + blk.log_singular_CO2_CO2_ratio[t, x], + blk.conc_CO2_equil_bulk[t, x], + blk.conc_CO2_bulk[t, x], + blk.log_conc_CO2_bulk[t, x], + ] + entangled_eqns = [ + blk.log_conc_interface_MEA_eqn[t, x], + blk.enhancement_factor_eqn1[t, x], + blk.log_singular_CO2_CO2_ratio_eqn[t, x], + blk.enhancement_factor_eqn2[t, x], + blk.conc_CO2_equil_bulk_eqn[t, x], + blk.conc_CO2_bulk_eqn[t, x], + blk.log_conc_CO2_bulk_eqn[t, x], + ] + long_var_list += entangled_vars + long_eqn_list += entangled_eqns + + tmp_blk = create_subsystem_block(long_eqn_list, long_var_list) + _sub_problem_scaling_suffix(blk, tmp_blk) + flags = _fix_vars([var for var in tmp_blk.input_vars.values()]) + assert degrees_of_freedom(tmp_blk) == 0 + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(tmp_blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + # if not check_optimal_termination(res): + # # IPOPT's solve failed, so let's try a sketchy iterative solution instead + # for t in blk.flowsheet().time: + # for x in blk.liquid_phase.length_domain: + # if x == blk.liquid_phase.length_domain.last(): + # continue + # zf = blk.liquid_phase.length_domain.next(x) + # for j in blk.log_diffus_liq_comp_list: + # blk.log_diffus_liq_comp[t, x, j].value = log(value(blk.liquid_phase.properties[t, x].diffus_phase_comp["Liq", j])) + # # CO2 is the only index + # blk.log_mass_transfer_coeff_liq[t, x, "CO2"].value = log(value(blk.mass_transfer_coeff_liq[t, x, "CO2"])) + + # calculate_variable_from_constraint( + # blk.log_hatta_number[t, x], blk.hatta_number_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + # ) + + + # if value(blk.conc_CO2_bulk[t, x]) < 1: + # blk.conc_interface_MEA[t, x].value = 0.95 + # blk.log_conc_interface_MEA[t, x].value = log(0.95) + # blk.log_singular_CO2_CO2_ratio[t, x].value = log(1.1) + # else: + # blk.conc_interface_MEA[t, x].value = 2 + # blk.log_conc_interface_MEA[t, x].value = log(2) + # blk.log_singular_CO2_CO2_ratio[t, x].value = log(0.9) + # blk.log_singular_MEA_CO2_ratio[t, x].value = value( + # log( + # (1 - blk.conc_interface_MEA[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_equil_bulk[t, x], + # blk.conc_CO2_equil_bulk_eqn[t, x], + # ) + # for k in range(8): + # # Solving Eq. 24 from Gaspar et. al. + # instant_E = value(exp(blk.log_instant_E_minus_one[t, x]) + 1) + # a = value(1 - instant_E) + # b = value(exp(blk.log_hatta_number[t, x]) * (blk.conc_CO2_equil_bulk[t, x] - 1)) + # c = value(instant_E - blk.conc_CO2_bulk[t, x]) + # Yplus = (-b + sqrt(b**2 - 4*a*c))/(2*a) + # Yminus = (-b - sqrt(b**2 - 4*a*c))/(2*a) + # if Yplus > 0: + # assert Yminus < 0 + # blk.conc_interface_MEA[t, x].value = Yplus**2 + # elif Yminus > 0: + # blk.conc_interface_MEA[t, x].value = Yminus**2 + # else: + # raise AssertionError + + # blk.log_conc_interface_MEA[t, x].value = log(blk.conc_interface_MEA[t, x].value) + + # # Use new value for conc_interface_MEA to calculate enhancement factor + # calculate_variable_from_constraint( + # blk.conc_CO2_equil_bulk[t, x], + # blk.conc_CO2_equil_bulk_eqn[t, x], + # ) + + # blk.log_singular_CO2_CO2_ratio[t, x].value = value( + # log( + # (1 - blk.conc_CO2_equil_bulk[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + # ) + + # enhancement_factor = value( + # exp( + # blk.log_hatta_number[t, x] + # + 0.5 * blk.log_conc_interface_MEA[t, x] + # + blk.log_singular_CO2_CO2_ratio[t, x] + # ) + # ) + # if enhancement_factor > 1: + # blk.log_enhancement_factor_minus_one[t, x].value = log(enhancement_factor - 1) + # else: + # blk.log_enhancement_factor_minus_one[t, x].value = -3 + # # Now update values with new value for enhancement factor + # calculate_variable_from_constraint( + # blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + # ) + # blk.log_singular_MEA_CO2_ratio[t, x].value = log( + # (1 - blk.conc_interface_MEA[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + + # # Resolve after (hopefully) better initialization + # assert degrees_of_freedom(tmp_blk) == 0 + # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + # res = solver_obj.solve(tmp_blk, tee=slc.tee, symbolic_solver_labels=True) + # assert_optimal_termination(res) + + + _restore_fixedness(flags) \ No newline at end of file diff --git a/idaes/models_extra/column_models/enhancement_factor_model_pseudo_second_order.py b/idaes/models_extra/column_models/enhancement_factor_model_pseudo_second_order.py new file mode 100644 index 0000000000..3ce4fa36ec --- /dev/null +++ b/idaes/models_extra/column_models/enhancement_factor_model_pseudo_second_order.py @@ -0,0 +1,695 @@ +# Import Python libraries +import logging +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd + +# Import Pyomo libraries +import pyomo.opt +from pyomo.environ import ( + Block, + ConcreteModel, + value, + Var, + Reals, + NonNegativeReals, + Param, + TransformationFactory, + Constraint, + Expression, + Objective, + SolverStatus, + TerminationCondition, + check_optimal_termination, + assert_optimal_termination, + exp, + log, + sqrt, + units as pyunits, + Set, + Reference +) +from pyomo.common.collections import ComponentSet, ComponentMap + +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.common.config import ConfigValue, Bool + +# Import IDAES Libraries +from idaes.core.util.constants import Constants as CONST +from idaes.models_extra.column_models.solvent_column import PackedColumnData + +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.initialization import _fix_vars, _restore_fixedness +from idaes.core import declare_process_block_class, FlowsheetBlock, StateBlock +from idaes.core.util.exceptions import InitializationError +from idaes.core.solvers.get_solver import get_solver +import idaes.logger as idaeslog + +from idaes.core.solvers import use_idaes_solver_configuration_defaults +import idaes.core.util.scaling as iscale +from pyomo.util.subsystems import ( + create_subsystem_block, +) +from idaes.core.solvers.petsc import ( + _sub_problem_scaling_suffix, +) + +def make_enhancement_factor_model(blk, lunits, kinetics="Putta"): + """ + Enhancement factor based liquid phase mass transfer model. + """ + assert kinetics in {"Luo", "Putta"} + blk.log_rate_constant_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Logarithm of rate constant for MEA mechanism", + initialize = 0, + ) + @blk.Constraint(blk.flowsheet().time, blk.liquid_phase.length_domain) + def log_rate_constant_MEA_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + if kinetics == "Putta": + # Putta, Svendsen, Knuutila 2017 Eqn. 42 + reduced_activation_energy_MEA = pyunits.convert( + 4936.6 * pyunits.K, to_units=lunits("temperature") + ) + preexponential_factor_MEA = pyunits.convert( + 3.1732e3 * ((pyunits.m) ** 6 / (pyunits.mol ** 2 * pyunits.s)), + to_units=1 / (lunits("time") * lunits("density_mole")**2), + ) + elif kinetics == "Luo": + reduced_activation_energy_MEA = pyunits.convert( + 4742.0 * pyunits.K, to_units=lunits("temperature") + ) + preexponential_factor_MEA = pyunits.convert( + 2.003e4 * ((pyunits.m) ** 6 / (pyunits.mol ** 2 * pyunits.s)), + to_units=1 / (lunits("time") * lunits("density_mole")**2), + ) + else: + return AssertionError + + log_preexponential_factor_MEA = log(value(preexponential_factor_MEA)) + + return b.log_rate_constant_MEA[t, x] == ( + log_preexponential_factor_MEA + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true["Liq", "MEA"] + - reduced_activation_energy_MEA / b.liquid_phase.properties[t, x].temperature + ) + + + blk.log_rate_constant_H2O = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Logarithm of rate constant for H2O mechanism", + initialize = 0, + ) + @blk.Constraint(blk.flowsheet().time, blk.liquid_phase.length_domain) + def log_rate_constant_H2O_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + if kinetics == "Putta": + # Putta, Svendsen, Knuutila 2017 Eqn. 42 + reduced_activation_energy_H2O = pyunits.convert( + 3900 * pyunits.K, to_units=lunits("temperature") + ) + preexponential_factor_H2O = pyunits.convert( + 1.0882e2 * ((pyunits.m) ** 6 / (pyunits.mol ** 2 * pyunits.s)), + to_units=1 / (lunits("time") * lunits("density_mole")**2), + ) + + elif kinetics == "Luo": + reduced_activation_energy_H2O = pyunits.convert( + 3110 * pyunits.K, to_units=lunits("temperature") + ) + preexponential_factor_H2O = pyunits.convert( + 4.147 * ((pyunits.m) ** 6 / (pyunits.mol ** 2 * pyunits.s)), + to_units=1 / (lunits("time") * lunits("density_mole")**2), + ) + else: + return AssertionError + + log_preexponential_factor_H2O = log(value(preexponential_factor_H2O)) + return b.log_rate_constant_H2O[t, x] == ( + log_preexponential_factor_H2O + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true["Liq", "H2O"] + - reduced_activation_energy_H2O / b.liquid_phase.properties[t, x].temperature + ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Second order rate constant [m3/(mol.s)]", + ) + def log_rate_constant(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return log( + exp(b.log_rate_constant_MEA[t, x]) + + exp(b.log_rate_constant_H2O[t, x]) + ) + + blk.log_hatta_number = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=3, + doc="Hatta number", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Hatta number constraint", + ) + def hatta_number_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return b.log_hatta_number[t, x] == 0.5 * ( + b.log_rate_constant[t, x] + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + + b.log_diffus_liq_comp[t, x, "CO2"] + ) - b.log_mass_transfer_coeff_liq[t, x, "CO2"] + + blk.conc_CO2_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=0.95, + units=pyunits.dimensionless, + bounds=(0, None), + doc="""Dimensionless concentration of CO2, + Driving force term where + Absorption implies conc_CO2_bulk < 1 and + Desorption implies conc_CO2_bulk > 1 """, + ) + blk.log_conc_CO2_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=-0.25, + units=pyunits.dimensionless, + doc="Natural logarithm of conc_CO2_bulk", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Equation defining log_conc_CO2_bulk", + ) + def log_conc_CO2_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return exp(b.log_conc_CO2_bulk[t, x]) == b.conc_CO2_bulk[t, x] + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Dimensionless concentration of CO2""", + ) + def conc_CO2_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + zf = b.liquid_phase.length_domain.next(x) + Pressure = pyunits.convert( + b.vapor_phase.properties[t, zf].pressure, + to_units=lunits("pressure"), + ) + return b.log_conc_CO2_bulk[t, x] + log( + ( + b.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + * Pressure + / b.psi[t, zf] + + exp(b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ]) + ) / lunits("density_mole") + ) == b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + log( + b.liquid_phase.properties[t, x].henry["Liq", "CO2"] + / b.psi[t, zf] + + 1 + ) + # return ( + # b.conc_CO2_bulk[t, x] * ( + # b.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + # * Pressure + # / b.psi[t, zf] + # + b.liquid_phase.properties[t, x].conc_mol_phase_comp_true["Liq", "CO2"] + # ) + # == b.liquid_phase.properties[t, x].conc_mol_phase_comp_true["Liq", "CO2"] + # * ( + # b.liquid_phase.properties[t, x].henry["Liq", "CO2"] + # / b.psi[t, zf] + # + 1 + # ) + # ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Instantaneous Enhancement factor""", + ) + def log_instant_E_minus_one(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + zf = b.liquid_phase.length_domain.next(x) + return ( + b.log_diffus_liq_comp[t, x, "MEA"] + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + + b.log_conc_CO2_bulk[t, x] + - log(2) + - b.log_diffus_liq_comp[t, x, "CO2"] + - b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + ) + # ====================================================================== + # Enhancement factor model + # Reference: Jozsef Gaspar,Philip Loldrup Fosbol, (2015) + + blk.conc_interface_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + bounds=(0, None), + initialize=0.95, + units=pyunits.dimensionless, + doc="""Dimensionless concentration of MEA + at interface """, + ) + + blk.log_conc_interface_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + bounds=(None, 100), + initialize=-0.05, + units=pyunits.dimensionless, + doc="""Logarithm of conc_interface_MEA""", + ) + # ============================================================================= + # ------------------------ ORIGINAL ----------------------------------- + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Dimensionless concentration of MEACOO-", + ) + def conc_interface_MEACOO(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return 1 + ( + b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + ) * (1 - b.conc_interface_MEA[t, x]) / ( + 2 + * b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEACOO_-" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEACOO_-" + ] + ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Dimensionless concentration of MEA+", + ) + def conc_interface_MEAH(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return 1 + ( + b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + ) * (1 - b.conc_interface_MEA[t, x]) / ( + 2 + * b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA_+" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA_+" + ] + ) + + # ============================================================================= + blk.conc_CO2_equil_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=1, + bounds=(0, None), + doc="""Dimensionless concentration of CO2 + at equilibrium with the bulk """, + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Constraint for dimensionless concentration of CO2 + at equilibrium with the bulk """, + ) + def conc_CO2_equil_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + b.conc_CO2_equil_bulk[t, x] * b.conc_interface_MEA[t, x] ** 2 + == b.conc_CO2_bulk[t, x] + * b.conc_interface_MEAH[t, x] + * b.conc_interface_MEACOO[t, x] + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Defining log_conc_interface_MEA", + ) + def log_conc_interface_MEA_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + exp(b.log_conc_interface_MEA[t, x]) + == b.conc_interface_MEA[t, x] + ) + + blk.log_singular_CO2_CO2_ratio = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=-3, + bounds=(-12, 12), + doc="""Logarithm of the ratio of one minus dimensionless concentration CO2 at equilibrium + with the bulk to one minus dimensionless concentration of CO2 actually in the bulk. + Individually, numerator and denominator go to zero at equilibrium, but the ratio + remains (hopefully) well defined mathematically, if not numerically.""", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Equation defining log_singular_CO2_CO2_ratio", + ) + def log_singular_CO2_CO2_ratio_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + # Constraint written presently to lose meaning when dimensionless concentrations approach 1 + # If this form doesn't work, we can try one with division in it + return exp(b.log_singular_CO2_CO2_ratio[t, x]) * ( + 1 - b.conc_CO2_bulk[t, x] + ) == (1 - b.conc_CO2_equil_bulk[t, x]) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Enhancement factor - function of Hatta number", + ) + def enhancement_factor_eqn1(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + b.log_enhancement_factor[t, x] + == b.log_hatta_number[t, x] + + 0.5 * b.log_conc_interface_MEA[t, x] + + b.log_singular_CO2_CO2_ratio[t, x] + ) + # return (b.enhancement_factor[t, x] * ( + # 1 - b.conc_CO2_bulk[t, x] + # ) == b.Hatta[t, x] * b.sqrt_conc_interface_MEA[t, x] * ( + # 1 - b.conc_CO2_equil_bulk[t, x] + # ) + + # blk.log_singular_MEA_CO2_ratio = Var( + # blk.flowsheet().time, + # blk.liquid_phase.length_domain, + # initialize=-5, + # bounds=( + # -12, + # 3, + # ), # Should be straight up nonpositive, but give it a bit of slack + # doc="""Logarithm of the ratio of one minus dimensionless concentration MEA at interface + # to one minus dimensionless concentration of CO2 at equilibrium with the bulk. + # Individually, numerator and denominator go to zero at equilibrium, but the ratio + # remains (hopefully) well defined mathematically, if not numerically.""", + # ) + + # @blk.Constraint( + # blk.flowsheet().time, + # blk.liquid_phase.length_domain, + # doc="Equation defining log_singular_MEA_CO2_ratio", + # ) + # def log_singular_MEA_CO2_ratio_eqn(b, t, x): + # if x == b.liquid_phase.length_domain.last(): + # return Constraint.Skip + # else: + # # Constraint written presently to lose meaning when dimensionless concentrations approach 1 + # # If this form doesn't work, we can try one with division in it + # return exp(b.log_singular_MEA_CO2_ratio[t, x]) * ( + # 1 - b.conc_CO2_bulk[t, x] + # ) == (1 - b.conc_interface_MEA[t, x]) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Enhancement factor - function of instantaneous enhancement factor", + ) + def enhancement_factor_eqn2(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + # return b.log_enhancement_factor_minus_one[t, x] == ( + # b.log_instant_E_minus_one[t, x] + # + b.log_singular_MEA_CO2_ratio[t, x] + # ) + return (b.enhancement_factor[t, x] - 1) * ( + 1 - b.conc_CO2_bulk[t, x] + ) == exp(b.log_instant_E_minus_one[t, x]) * (1 - b.conc_interface_MEA[t, x]) + + enhancement_factor_vars = [ + blk.log_enhancement_factor, + blk.log_hatta_number, + blk.conc_CO2_bulk, + blk.log_conc_CO2_bulk, + blk.conc_interface_MEA, + blk.log_conc_interface_MEA, + blk.conc_CO2_equil_bulk, + blk.log_singular_CO2_CO2_ratio, + blk.log_rate_constant_MEA, + blk.log_rate_constant_H2O, + ] + + enhancement_factor_constraints = [ + blk.hatta_number_eqn, + blk.log_conc_CO2_bulk_eqn, + blk.conc_CO2_bulk_eqn, + blk.conc_CO2_equil_bulk_eqn, + blk.log_conc_interface_MEA_eqn, + blk.log_singular_CO2_CO2_ratio_eqn, + blk.enhancement_factor_eqn1, + blk.enhancement_factor_eqn2, + blk.log_rate_constant_MEA_eqn, + blk.log_rate_constant_H2O_eqn, + ] + + return enhancement_factor_vars, enhancement_factor_constraints + +def initialize_enhancement_factor_model( + blk, + state_args=None, + outlvl=idaeslog.NOTSET, + optarg=None, + solver=None, + ): + # Set up logger for initialization and solve + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") + # Set solver options + if optarg is None: + optarg = {} + + solver_obj = get_solver(solver, optarg) + + long_var_list = [] + long_eqn_list = [] + + for t in blk.flowsheet().time: + for x in blk.liquid_phase.length_domain: + if x == blk.liquid_phase.length_domain.last(): + continue + zf = blk.liquid_phase.length_domain.next(x) + + calculate_variable_from_constraint( + blk.log_rate_constant_MEA[t, x], blk.log_rate_constant_MEA_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_rate_constant_H2O[t, x], blk.log_rate_constant_H2O_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_hatta_number[t, x], blk.hatta_number_eqn[t, x] + ) + + if value(blk.conc_CO2_bulk[t, x]) < 1: + blk.conc_interface_MEA[t, x].value = 0.95 + blk.log_conc_interface_MEA[t, x].value = log(0.95) + blk.log_singular_CO2_CO2_ratio[t, x].value = log(1.1) + else: + blk.conc_interface_MEA[t, x].value = 1.05 + blk.log_conc_interface_MEA[t, x].value = log(1.05) + blk.log_singular_CO2_CO2_ratio[t, x].value = log(0.9) + calculate_variable_from_constraint( + blk.conc_CO2_equil_bulk[t, x], + blk.conc_CO2_equil_bulk_eqn[t, x], + ) + + entangled_vars = [ + blk.conc_interface_MEA[t, x], + blk.log_conc_interface_MEA[t, x], + blk.log_enhancement_factor[t, x], + blk.log_singular_CO2_CO2_ratio[t, x], + blk.conc_CO2_equil_bulk[t, x], + blk.conc_CO2_bulk[t, x], + blk.log_conc_CO2_bulk[t, x], + ] + entangled_eqns = [ + blk.log_conc_interface_MEA_eqn[t, x], + blk.enhancement_factor_eqn1[t, x], + blk.log_singular_CO2_CO2_ratio_eqn[t, x], + blk.enhancement_factor_eqn2[t, x], + blk.conc_CO2_equil_bulk_eqn[t, x], + blk.conc_CO2_bulk_eqn[t, x], + blk.log_conc_CO2_bulk_eqn[t, x], + ] + long_var_list += entangled_vars + long_eqn_list += entangled_eqns + + tmp_blk = create_subsystem_block(long_eqn_list, long_var_list) + _sub_problem_scaling_suffix(blk, tmp_blk) + flags = _fix_vars([var for var in tmp_blk.input_vars.values()]) + assert degrees_of_freedom(tmp_blk) == 0 + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(tmp_blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + # if not check_optimal_termination(res): + # # IPOPT's solve failed, so let's try a sketchy iterative solution instead + # for t in blk.flowsheet().time: + # for x in blk.liquid_phase.length_domain: + # if x == blk.liquid_phase.length_domain.last(): + # continue + # zf = blk.liquid_phase.length_domain.next(x) + # for j in blk.log_diffus_liq_comp_list: + # blk.log_diffus_liq_comp[t, x, j].value = log(value(blk.liquid_phase.properties[t, x].diffus_phase_comp["Liq", j])) + # # CO2 is the only index + # blk.log_mass_transfer_coeff_liq[t, x, "CO2"].value = log(value(blk.mass_transfer_coeff_liq[t, x, "CO2"])) + + # calculate_variable_from_constraint( + # blk.log_hatta_number[t, x], blk.hatta_number_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + # ) + + + # if value(blk.conc_CO2_bulk[t, x]) < 1: + # blk.conc_interface_MEA[t, x].value = 0.95 + # blk.log_conc_interface_MEA[t, x].value = log(0.95) + # blk.log_singular_CO2_CO2_ratio[t, x].value = log(1.1) + # else: + # blk.conc_interface_MEA[t, x].value = 2 + # blk.log_conc_interface_MEA[t, x].value = log(2) + # blk.log_singular_CO2_CO2_ratio[t, x].value = log(0.9) + # blk.log_singular_MEA_CO2_ratio[t, x].value = value( + # log( + # (1 - blk.conc_interface_MEA[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_equil_bulk[t, x], + # blk.conc_CO2_equil_bulk_eqn[t, x], + # ) + # for k in range(8): + # # Solving Eq. 24 from Gaspar et. al. + # instant_E = value(exp(blk.log_instant_E_minus_one[t, x]) + 1) + # a = value(1 - instant_E) + # b = value(exp(blk.log_hatta_number[t, x]) * (blk.conc_CO2_equil_bulk[t, x] - 1)) + # c = value(instant_E - blk.conc_CO2_bulk[t, x]) + # Yplus = (-b + sqrt(b**2 - 4*a*c))/(2*a) + # Yminus = (-b - sqrt(b**2 - 4*a*c))/(2*a) + # if Yplus > 0: + # assert Yminus < 0 + # blk.conc_interface_MEA[t, x].value = Yplus**2 + # elif Yminus > 0: + # blk.conc_interface_MEA[t, x].value = Yminus**2 + # else: + # raise AssertionError + + # blk.log_conc_interface_MEA[t, x].value = log(blk.conc_interface_MEA[t, x].value) + + # # Use new value for conc_interface_MEA to calculate enhancement factor + # calculate_variable_from_constraint( + # blk.conc_CO2_equil_bulk[t, x], + # blk.conc_CO2_equil_bulk_eqn[t, x], + # ) + + # blk.log_singular_CO2_CO2_ratio[t, x].value = value( + # log( + # (1 - blk.conc_CO2_equil_bulk[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + # ) + + # enhancement_factor = value( + # exp( + # blk.log_hatta_number[t, x] + # + 0.5 * blk.log_conc_interface_MEA[t, x] + # + blk.log_singular_CO2_CO2_ratio[t, x] + # ) + # ) + # if enhancement_factor > 1: + # blk.log_enhancement_factor_minus_one[t, x].value = log(enhancement_factor - 1) + # else: + # blk.log_enhancement_factor_minus_one[t, x].value = -3 + # # Now update values with new value for enhancement factor + # calculate_variable_from_constraint( + # blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + # ) + # blk.log_singular_MEA_CO2_ratio[t, x].value = log( + # (1 - blk.conc_interface_MEA[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + + # # Resolve after (hopefully) better initialization + # assert degrees_of_freedom(tmp_blk) == 0 + # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + # res = solver_obj.solve(tmp_blk, tee=slc.tee, symbolic_solver_labels=True) + # assert_optimal_termination(res) + + + _restore_fixedness(flags) \ No newline at end of file diff --git a/idaes/models_extra/column_models/enhancement_factor_model_third_order.py b/idaes/models_extra/column_models/enhancement_factor_model_third_order.py new file mode 100644 index 0000000000..7649edeec6 --- /dev/null +++ b/idaes/models_extra/column_models/enhancement_factor_model_third_order.py @@ -0,0 +1,622 @@ +# Import Python libraries +import logging +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd + +# Import Pyomo libraries +import pyomo.opt +from pyomo.environ import ( + Block, + ConcreteModel, + value, + Var, + Reals, + NonNegativeReals, + Param, + TransformationFactory, + Constraint, + Expression, + Objective, + SolverStatus, + TerminationCondition, + check_optimal_termination, + assert_optimal_termination, + exp, + log, + sqrt, + units as pyunits, + Set, + Reference +) +from pyomo.common.collections import ComponentSet, ComponentMap + +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.common.config import ConfigValue, Bool + +# Import IDAES Libraries +from idaes.core.util.constants import Constants as CONST +from idaes.models_extra.column_models.solvent_column import PackedColumnData + +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.initialization import _fix_vars, _restore_fixedness +from idaes.core import declare_process_block_class, FlowsheetBlock, StateBlock +from idaes.core.util.exceptions import InitializationError +from idaes.core.solvers.get_solver import get_solver +import idaes.logger as idaeslog + +from idaes.core.solvers import use_idaes_solver_configuration_defaults +import idaes.core.util.scaling as iscale +from pyomo.util.subsystems import ( + create_subsystem_block, +) +from idaes.core.solvers.petsc import ( + _sub_problem_scaling_suffix, +) + +def make_enhancement_factor_model(blk, lunits): + """ + Enhancement factor based liquid phase mass transfer model. + """ + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Second order rate constant [m3/(mol.s)]", + ) + def log_rate_constant(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + T = pyunits.convert( + b.liquid_phase.properties[t, x].temperature, to_units=pyunits.K + ) + # C_MEA = pyunits.convert( + # b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + # "Liq", "MEA" + # ], + # to_units=pyunits.mol / pyunits.m ** 3, + # ) + + # Putta, Svendsen, Knuutila 2017 Eqn. 42 + return log( + pyunits.convert( + 3.1732e9 + * exp(-4936.6 * pyunits.K / T) + * ((pyunits.m) ** 6 / (pyunits.kmol ** 2 * pyunits.s)), + to_units=1 / (lunits("time") * lunits("density_mole") ** 2), + ) + * lunits("time") + * lunits("density_mole") ** 2 + ) + + blk.log_hatta_number = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=3, + doc="Hatta number", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Hatta number constraint", + ) + def hatta_number_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return b.log_hatta_number[t, x] == 0.5 * ( + b.log_rate_constant[t, x] + + 2 * b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + + b.log_diffus_liq_comp[t, x, "CO2"] + ) - b.log_mass_transfer_coeff_liq[t, x, "CO2"] + + blk.conc_CO2_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=0.95, + units=pyunits.dimensionless, + bounds=(0, None), + doc="""Dimensionless concentration of CO2, + Driving force term where + Absorption implies conc_CO2_bulk < 1 and + Desorption implies conc_CO2_bulk > 1 """, + ) + blk.log_conc_CO2_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=-0.25, + units=pyunits.dimensionless, + doc="Natural logarithm of conc_CO2_bulk", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Equation defining log_conc_CO2_bulk", + ) + def log_conc_CO2_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return exp(b.log_conc_CO2_bulk[t, x]) == b.conc_CO2_bulk[t, x] + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Dimensionless concentration of CO2""", + ) + def conc_CO2_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + zf = b.liquid_phase.length_domain.next(x) + Pressure = pyunits.convert( + b.vapor_phase.properties[t, zf].pressure, + to_units=lunits("pressure"), + ) + return b.log_conc_CO2_bulk[t, x] + log( + ( + b.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + * Pressure + / b.psi[t, zf] + + b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + ) / lunits("density_mole") + ) == b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + log( + b.liquid_phase.properties[t, x].henry["Liq", "CO2"] + / b.psi[t, zf] + + 1 + ) + # return ( + # b.conc_CO2_bulk[t, x] * ( + # b.vapor_phase.properties[t, zf].mole_frac_comp["CO2"] + # * Pressure + # / b.psi[t, zf] + # + b.liquid_phase.properties[t, x].conc_mol_phase_comp_true["Liq", "CO2"] + # ) + # == b.liquid_phase.properties[t, x].conc_mol_phase_comp_true["Liq", "CO2"] + # * ( + # b.liquid_phase.properties[t, x].henry["Liq", "CO2"] + # / b.psi[t, zf] + # + 1 + # ) + # ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Instantaneous Enhancement factor""", + ) + def log_instant_E_minus_one(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + zf = b.liquid_phase.length_domain.next(x) + return ( + b.log_diffus_liq_comp[t, x, "MEA"] + + b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + + b.log_conc_CO2_bulk[t, x] + - log(2) + - b.log_diffus_liq_comp[t, x, "CO2"] + - b.liquid_phase.properties[t, x].log_conc_mol_phase_comp_true[ + "Liq", "CO2" + ] + ) + # ====================================================================== + # Enhancement factor model + # Reference: Jozsef Gaspar,Philip Loldrup Fosbol, (2015) + + blk.conc_interface_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + bounds=(0, None), + initialize=0.95, + units=pyunits.dimensionless, + doc="""Dimensionless concentration of MEA + at interface """, + ) + + blk.log_conc_interface_MEA = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + bounds=(None, 100), + initialize=-0.05, + units=pyunits.dimensionless, + doc="""Logarithm of conc_interface_MEA""", + ) + # ============================================================================= + # ------------------------ ORIGINAL ----------------------------------- + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Dimensionless concentration of MEACOO-", + ) + def conc_interface_MEACOO(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return 1 + ( + b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + ) * (1 - b.conc_interface_MEA[t, x]) / ( + 2 + * b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEACOO_-" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEACOO_-" + ] + ) + + @blk.Expression( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Dimensionless concentration of MEA+", + ) + def conc_interface_MEAH(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Expression.Skip + else: + return 1 + ( + b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA" + ] + ) * (1 - b.conc_interface_MEA[t, x]) / ( + 2 + * b.liquid_phase.properties[t, x].diffus_phase_comp_true[ + "Liq", "MEA_+" + ] + * b.liquid_phase.properties[t, x].conc_mol_phase_comp_true[ + "Liq", "MEA_+" + ] + ) + + # ============================================================================= + blk.conc_CO2_equil_bulk = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=1, + bounds=(0, None), + doc="""Dimensionless concentration of CO2 + at equilibrium with the bulk """, + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="""Constraint for dimensionless concentration of CO2 + at equilibrium with the bulk """, + ) + def conc_CO2_equil_bulk_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + b.conc_CO2_equil_bulk[t, x] * b.conc_interface_MEA[t, x] ** 2 + == b.conc_CO2_bulk[t, x] + * b.conc_interface_MEAH[t, x] + * b.conc_interface_MEACOO[t, x] + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Defining log_conc_interface_MEA", + ) + def log_conc_interface_MEA_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + exp(b.log_conc_interface_MEA[t, x]) + == b.conc_interface_MEA[t, x] + ) + + blk.log_singular_CO2_CO2_ratio = Var( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + initialize=-3, + bounds=(-12, 12), + doc="""Logarithm of the ratio of one minus dimensionless concentration CO2 at equilibrium + with the bulk to one minus dimensionless concentration of CO2 actually in the bulk. + Individually, numerator and denominator go to zero at equilibrium, but the ratio + remains (hopefully) well defined mathematically, if not numerically.""", + ) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Equation defining log_singular_CO2_CO2_ratio", + ) + def log_singular_CO2_CO2_ratio_eqn(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + # Constraint written presently to lose meaning when dimensionless concentrations approach 1 + # If this form doesn't work, we can try one with division in it + return exp(b.log_singular_CO2_CO2_ratio[t, x]) * ( + 1 - b.conc_CO2_bulk[t, x] + ) == (1 - b.conc_CO2_equil_bulk[t, x]) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Enhancement factor - function of Hatta number", + ) + def enhancement_factor_eqn1(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + return ( + b.log_enhancement_factor[t, x] + == b.log_hatta_number[t, x] + + b.log_conc_interface_MEA[t, x] + + b.log_singular_CO2_CO2_ratio[t, x] + ) + # return (b.enhancement_factor[t, x] * ( + # 1 - b.conc_CO2_bulk[t, x] + # ) == b.Hatta[t, x] * b.sqrt_conc_interface_MEA[t, x] * ( + # 1 - b.conc_CO2_equil_bulk[t, x] + # ) + + # blk.log_singular_MEA_CO2_ratio = Var( + # blk.flowsheet().time, + # blk.liquid_phase.length_domain, + # initialize=-5, + # bounds=( + # -12, + # 3, + # ), # Should be straight up nonpositive, but give it a bit of slack + # doc="""Logarithm of the ratio of one minus dimensionless concentration MEA at interface + # to one minus dimensionless concentration of CO2 at equilibrium with the bulk. + # Individually, numerator and denominator go to zero at equilibrium, but the ratio + # remains (hopefully) well defined mathematically, if not numerically.""", + # ) + + # @blk.Constraint( + # blk.flowsheet().time, + # blk.liquid_phase.length_domain, + # doc="Equation defining log_singular_MEA_CO2_ratio", + # ) + # def log_singular_MEA_CO2_ratio_eqn(b, t, x): + # if x == b.liquid_phase.length_domain.last(): + # return Constraint.Skip + # else: + # # Constraint written presently to lose meaning when dimensionless concentrations approach 1 + # # If this form doesn't work, we can try one with division in it + # return exp(b.log_singular_MEA_CO2_ratio[t, x]) * ( + # 1 - b.conc_CO2_bulk[t, x] + # ) == (1 - b.conc_interface_MEA[t, x]) + + @blk.Constraint( + blk.flowsheet().time, + blk.liquid_phase.length_domain, + doc="Enhancement factor - function of instantaneous enhancement factor", + ) + def enhancement_factor_eqn2(b, t, x): + if x == b.liquid_phase.length_domain.last(): + return Constraint.Skip + else: + # return b.log_enhancement_factor_minus_one[t, x] == ( + # b.log_instant_E_minus_one[t, x] + # + b.log_singular_MEA_CO2_ratio[t, x] + # ) + return (b.enhancement_factor[t, x] - 1) * ( + 1 - b.conc_CO2_bulk[t, x] + ) == exp(b.log_instant_E_minus_one[t, x]) * (1 - b.conc_interface_MEA[t, x]) + + enhancement_factor_vars = [ + blk.log_enhancement_factor, + blk.log_hatta_number, + blk.conc_CO2_bulk, + blk.log_conc_CO2_bulk, + blk.conc_interface_MEA, + blk.log_conc_interface_MEA, + blk.conc_CO2_equil_bulk, + blk.log_singular_CO2_CO2_ratio, + ] + + enhancement_factor_constraints = [ + blk.hatta_number_eqn, + blk.log_conc_CO2_bulk_eqn, + blk.conc_CO2_bulk_eqn, + blk.conc_CO2_equil_bulk_eqn, + blk.log_conc_interface_MEA_eqn, + blk.log_singular_CO2_CO2_ratio_eqn, + blk.enhancement_factor_eqn1, + blk.enhancement_factor_eqn2, + ] + + return enhancement_factor_vars, enhancement_factor_constraints + +def initialize_enhancement_factor_model( + blk, + state_args=None, + outlvl=idaeslog.NOTSET, + optarg=None, + solver=None, + ): + # Set up logger for initialization and solve + init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") + # Set solver options + if optarg is None: + optarg = {} + + solver_obj = get_solver(solver, optarg) + + long_var_list = [] + long_eqn_list = [] + + for t in blk.flowsheet().time: + for x in blk.liquid_phase.length_domain: + if x == blk.liquid_phase.length_domain.last(): + continue + zf = blk.liquid_phase.length_domain.next(x) + calculate_variable_from_constraint( + blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + ) + calculate_variable_from_constraint( + blk.log_hatta_number[t, x], blk.hatta_number_eqn[t, x] + ) + + if value(blk.conc_CO2_bulk[t, x]) < 1: + blk.conc_interface_MEA[t, x].value = 0.95 + blk.log_conc_interface_MEA[t, x].value = log(0.95) + blk.log_singular_CO2_CO2_ratio[t, x].value = log(1.1) + else: + blk.conc_interface_MEA[t, x].value = 1.05 + blk.log_conc_interface_MEA[t, x].value = log(1.05) + blk.log_singular_CO2_CO2_ratio[t, x].value = log(0.9) + calculate_variable_from_constraint( + blk.conc_CO2_equil_bulk[t, x], + blk.conc_CO2_equil_bulk_eqn[t, x], + ) + + entangled_vars = [ + blk.conc_interface_MEA[t, x], + blk.log_conc_interface_MEA[t, x], + blk.log_enhancement_factor[t, x], + blk.log_singular_CO2_CO2_ratio[t, x], + blk.conc_CO2_equil_bulk[t, x], + blk.conc_CO2_bulk[t, x], + blk.log_conc_CO2_bulk[t, x], + ] + entangled_eqns = [ + blk.log_conc_interface_MEA_eqn[t, x], + blk.enhancement_factor_eqn1[t, x], + blk.log_singular_CO2_CO2_ratio_eqn[t, x], + blk.enhancement_factor_eqn2[t, x], + blk.conc_CO2_equil_bulk_eqn[t, x], + blk.conc_CO2_bulk_eqn[t, x], + blk.log_conc_CO2_bulk_eqn[t, x], + ] + long_var_list += entangled_vars + long_eqn_list += entangled_eqns + + tmp_blk = create_subsystem_block(long_eqn_list, long_var_list) + _sub_problem_scaling_suffix(blk, tmp_blk) + flags = _fix_vars([var for var in tmp_blk.input_vars.values()]) + assert degrees_of_freedom(tmp_blk) == 0 + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(tmp_blk, tee=slc.tee, symbolic_solver_labels=True) + assert_optimal_termination(res) + # if not check_optimal_termination(res): + # # IPOPT's solve failed, so let's try a sketchy iterative solution instead + # for t in blk.flowsheet().time: + # for x in blk.liquid_phase.length_domain: + # if x == blk.liquid_phase.length_domain.last(): + # continue + # zf = blk.liquid_phase.length_domain.next(x) + # for j in blk.log_diffus_liq_comp_list: + # blk.log_diffus_liq_comp[t, x, j].value = log(value(blk.liquid_phase.properties[t, x].diffus_phase_comp["Liq", j])) + # # CO2 is the only index + # blk.log_mass_transfer_coeff_liq[t, x, "CO2"].value = log(value(blk.mass_transfer_coeff_liq[t, x, "CO2"])) + + # calculate_variable_from_constraint( + # blk.log_hatta_number[t, x], blk.hatta_number_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + # ) + + + # if value(blk.conc_CO2_bulk[t, x]) < 1: + # blk.conc_interface_MEA[t, x].value = 0.95 + # blk.log_conc_interface_MEA[t, x].value = log(0.95) + # blk.log_singular_CO2_CO2_ratio[t, x].value = log(1.1) + # else: + # blk.conc_interface_MEA[t, x].value = 2 + # blk.log_conc_interface_MEA[t, x].value = log(2) + # blk.log_singular_CO2_CO2_ratio[t, x].value = log(0.9) + # blk.log_singular_MEA_CO2_ratio[t, x].value = value( + # log( + # (1 - blk.conc_interface_MEA[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_equil_bulk[t, x], + # blk.conc_CO2_equil_bulk_eqn[t, x], + # ) + # for k in range(8): + # # Solving Eq. 24 from Gaspar et. al. + # instant_E = value(exp(blk.log_instant_E_minus_one[t, x]) + 1) + # a = value(1 - instant_E) + # b = value(exp(blk.log_hatta_number[t, x]) * (blk.conc_CO2_equil_bulk[t, x] - 1)) + # c = value(instant_E - blk.conc_CO2_bulk[t, x]) + # Yplus = (-b + sqrt(b**2 - 4*a*c))/(2*a) + # Yminus = (-b - sqrt(b**2 - 4*a*c))/(2*a) + # if Yplus > 0: + # assert Yminus < 0 + # blk.conc_interface_MEA[t, x].value = Yplus**2 + # elif Yminus > 0: + # blk.conc_interface_MEA[t, x].value = Yminus**2 + # else: + # raise AssertionError + + # blk.log_conc_interface_MEA[t, x].value = log(blk.conc_interface_MEA[t, x].value) + + # # Use new value for conc_interface_MEA to calculate enhancement factor + # calculate_variable_from_constraint( + # blk.conc_CO2_equil_bulk[t, x], + # blk.conc_CO2_equil_bulk_eqn[t, x], + # ) + + # blk.log_singular_CO2_CO2_ratio[t, x].value = value( + # log( + # (1 - blk.conc_CO2_equil_bulk[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + # ) + + # enhancement_factor = value( + # exp( + # blk.log_hatta_number[t, x] + # + 0.5 * blk.log_conc_interface_MEA[t, x] + # + blk.log_singular_CO2_CO2_ratio[t, x] + # ) + # ) + # if enhancement_factor > 1: + # blk.log_enhancement_factor_minus_one[t, x].value = log(enhancement_factor - 1) + # else: + # blk.log_enhancement_factor_minus_one[t, x].value = -3 + # # Now update values with new value for enhancement factor + # calculate_variable_from_constraint( + # blk.log_conc_CO2_bulk[t, x], blk.conc_CO2_bulk_eqn[t, x] + # ) + # calculate_variable_from_constraint( + # blk.conc_CO2_bulk[t, x], blk.log_conc_CO2_bulk_eqn[t, x] + # ) + # blk.log_singular_MEA_CO2_ratio[t, x].value = log( + # (1 - blk.conc_interface_MEA[t, x]) / (1 - blk.conc_CO2_bulk[t, x]) + # ) + + # # Resolve after (hopefully) better initialization + # assert degrees_of_freedom(tmp_blk) == 0 + # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + # res = solver_obj.solve(tmp_blk, tee=slc.tee, symbolic_solver_labels=True) + # assert_optimal_termination(res) + + + _restore_fixedness(flags) \ No newline at end of file diff --git a/idaes/models_extra/column_models/properties/MEA_solvent.py b/idaes/models_extra/column_models/properties/MEA_solvent.py index 264b6373ef..25e8d9cd48 100644 --- a/idaes/models_extra/column_models/properties/MEA_solvent.py +++ b/idaes/models_extra/column_models/properties/MEA_solvent.py @@ -255,7 +255,7 @@ def return_expression(b, p, j, T): class PressureSatSolvent: - # Method for calculating saturation pressure ofsolvents + # Method for calculating saturation pressure of solvents @staticmethod def build_parameters(cobj): cobj.pressure_sat_comp_coeff_1 = Var( @@ -331,17 +331,27 @@ def return_expression(b, cobj, T): return pyunits.convert(vol_mol, units.VOLUME / units.AMOUNT) - -class VolMolCO2: - # Weiland Method for calculating molar volume of dissolved CO2 [2] - +class VolMolMEA: + # Molar volume of MEA as calculated by Morgan [2] @staticmethod def build_parameters(cobj): - cobj.vol_mol_liq_comp_coeff_a = Var( - doc="Parameter a for liquid phase molar volume", - units=pyunits.mL / pyunits.mol, + cobj.dens_mol_liq_comp_coeff_1 = Var( + doc="Parameter 1 for liquid phase molar density", + units=pyunits.g / pyunits.mL / pyunits.K**2, ) - set_param_from_config(cobj, param="vol_mol_liq_comp_coeff", index="a") + set_param_from_config(cobj, param="dens_mol_liq_comp_coeff", index="1") + + cobj.dens_mol_liq_comp_coeff_2 = Var( + doc="Parameter 2 for liquid phase molar density", + units=pyunits.g / pyunits.mL / pyunits.K, + ) + set_param_from_config(cobj, param="dens_mol_liq_comp_coeff", index="2") + + cobj.dens_mol_liq_comp_coeff_3 = Var( + doc="Parameter 3 for liquid phase molar density", + units=pyunits.g / pyunits.mL, + ) + set_param_from_config(cobj, param="dens_mol_liq_comp_coeff", index="3") cobj.vol_mol_liq_comp_coeff_b = Var( doc="Parameter b for liquid phase molar volume", @@ -355,6 +365,38 @@ def build_parameters(cobj): ) set_param_from_config(cobj, param="vol_mol_liq_comp_coeff", index="c") + @staticmethod + def return_expression(b, cobj, T): + units = b.params.get_metadata().derived_units + T = pyunits.convert(T, to_units=pyunits.K) + x = b.mole_frac_comp + + rho = ( + cobj.dens_mol_liq_comp_coeff_1 * T**2 + + cobj.dens_mol_liq_comp_coeff_2 * T + + cobj.dens_mol_liq_comp_coeff_3 + ) + vol_mol_pure = pyunits.convert(cobj.mw / rho, units.VOLUME / units.AMOUNT) + + vol_mol_interaction = x["H2O"] * ( + cobj.vol_mol_liq_comp_coeff_b + + cobj.vol_mol_liq_comp_coeff_c * x["MEA"] + ) + + + return vol_mol_pure + pyunits.convert(vol_mol_interaction, units.VOLUME / units.AMOUNT) + +class VolMolCO2: + # Weiland Method for calculating molar volume of dissolved CO2 [2] + + @staticmethod + def build_parameters(cobj): + cobj.vol_mol_liq_comp_coeff_a = Var( + doc="Parameter a for liquid phase molar volume", + units=pyunits.mL / pyunits.mol, + ) + set_param_from_config(cobj, param="vol_mol_liq_comp_coeff", index="a") + cobj.vol_mol_liq_comp_coeff_d = Var( doc="Parameter d for liquid phase molar volume", units=pyunits.mL / pyunits.mol, @@ -373,10 +415,6 @@ def return_expression(b, cobj, T): vol_mol = ( cobj.vol_mol_liq_comp_coeff_a - + (cobj.vol_mol_liq_comp_coeff_b + cobj.vol_mol_liq_comp_coeff_c * x["MEA"]) - * x["MEA"] - * x["H2O"] - / x["CO2"] + (cobj.vol_mol_liq_comp_coeff_d + cobj.vol_mol_liq_comp_coeff_e * x["MEA"]) * x["MEA"] ) @@ -904,7 +942,7 @@ def calculate_scaling_factors(b, rblock): "diffus_phase_comp": {"Liq": DiffusMEA}, "enth_mol_liq_comp": EnthMolSolvent, "pressure_sat_comp": PressureSatSolvent, - "vol_mol_liq_comp": VolMolSolvent, + "vol_mol_liq_comp": VolMolMEA, "parameter_data": { "mw": (0.06108, pyunits.kg / pyunits.mol), "cp_mass_liq_comp_coeff": { @@ -919,6 +957,10 @@ def calculate_scaling_factors(b, rblock): "2": (-4.51417e-4, pyunits.g / pyunits.mL / pyunits.K), "3": (1.19451, pyunits.g / pyunits.mL), }, + "vol_mol_liq_comp_coeff": { + "b": (-2.2642, pyunits.mL / pyunits.mol), # [2] + "c": (3.0059, pyunits.mL / pyunits.mol), + }, "dh_vap": 58000, # [3] "diffus_phase_comp_coeff": { "1": -13.275, @@ -965,8 +1007,6 @@ def calculate_scaling_factors(b, rblock): }, "vol_mol_liq_comp_coeff": { "a": (10.2074, pyunits.mL / pyunits.mol), # [2] - "b": (-2.2642, pyunits.mL / pyunits.mol), - "c": (3.0059, pyunits.mL / pyunits.mol), "d": (207, pyunits.mL / pyunits.mol), "e": (-563.3701, pyunits.mL / pyunits.mol), }, diff --git a/idaes/models_extra/column_models/solvent_column.py b/idaes/models_extra/column_models/solvent_column.py index e13c2111b5..4ef8408995 100644 --- a/idaes/models_extra/column_models/solvent_column.py +++ b/idaes/models_extra/column_models/solvent_column.py @@ -372,12 +372,14 @@ def column_cross_section_area_eqn(blk): expr=4 * self.eps_ref / self.packing_specific_area, doc="Hydraulic diameter" ) + # TODO Change this to specific_interfacial_area self.area_interfacial = Var( self.flowsheet().time, self.vapor_phase.length_domain, initialize=0.9, - units=(pyunits.m) ** 2 / (pyunits.m) ** 3, - doc="Specific interfacial area", + units=lunits("length") ** 2 / lunits("length") ** 3, + # TODO this description is wrong + doc="Interface area between vapor and liquid per unit of column volume", ) # Liquid and vapor holdups @@ -437,7 +439,7 @@ def liquid_phase_area(blk, t, x): @self.Constraint( self.flowsheet().time, self.liquid_phase.length_domain, - doc="Mechanical equilibrium constraint", + doc="Mechanical equilibruim constraint", ) def mechanical_equilibrium(blk, t, x): if x == self.liquid_phase.length_domain.last(): @@ -454,6 +456,7 @@ def mechanical_equilibrium(blk, t, x): self.flowsheet().time, self.vapor_phase.length_domain, equilibrium_comp, + bounds=(0, None), units=( lunits("amount") / lunits("pressure") @@ -463,18 +466,29 @@ def mechanical_equilibrium(blk, t, x): doc="Vapor phase mass transfer coefficient", ) - # Equilibrium partial pressure of components at interface - self.pressure_equil = Var( + # Equilibruim partial pressure of components at interface + # self.pressure_equil = Var( + # self.flowsheet().time, + # self.vapor_phase.length_domain, + # equilibrium_comp, + # domain=NonNegativeReals, + # initialize=500, + # units=lunits("pressure"), + # doc="Equilibruim pressure of components at interface", + # ) + + # Mass transfer constraints + # "mass_transfer" is a bad name for these variables since they're written with a mole basis. "material" is better + self.mass_transfer_driving_force = Var( self.flowsheet().time, self.vapor_phase.length_domain, equilibrium_comp, - domain=NonNegativeReals, - initialize=500, + initialize=0, + domain=Reals, units=lunits("pressure"), - doc="Equilibrium pressure of components at interface", + doc="Generalized driving force for mass transfer" ) - # Mass transfer constraints self.interphase_mass_transfer = Var( self.flowsheet().time, self.liquid_phase.length_domain, @@ -501,14 +515,7 @@ def interphase_mass_transfer_eqn(blk, t, x, j): blk.mass_transfer_coeff_vap[t, x, j] * blk.area_interfacial[t, x] * blk.area_column - * ( - blk.vapor_phase.properties[t, x].mole_frac_comp[j] - * pyunits.convert( - blk.vapor_phase.properties[t, x].pressure, - to_units=lunits("pressure"), - ) - - blk.pressure_equil[t, x, j] - ) + * blk.mass_transfer_driving_force[t, x, j] ) else: return blk.interphase_mass_transfer[t, x, j] == 0.0 @@ -531,6 +538,9 @@ def liquid_mass_transfer_eqn(blk, t, x, j): else: return blk.liquid_phase.mass_transfer_term[t, x, "Liq", j] == 0.0 + vunits = ( + self.config.vapor_phase.property_package.get_metadata().get_derived_units + ) @self.Constraint( self.flowsheet().time, self.vapor_phase.length_domain, @@ -543,10 +553,10 @@ def vapor_mass_transfer_eqn(blk, t, x, j): elif j in equilibrium_comp: return ( pyunits.convert( - blk.vapor_phase.mass_transfer_term[t, x, "Vap", j], - to_units=lunits("amount") / lunits("time") / lunits("length"), + -blk.interphase_mass_transfer[t, x, j], + to_units=vunits("amount") / vunits("time") / vunits("length"), ) - == -blk.interphase_mass_transfer[t, x, j] + == blk.vapor_phase.mass_transfer_term[t, x, "Vap", j] ) else: return blk.vapor_phase.mass_transfer_term[t, x, "Vap", j] == 0.0 @@ -558,7 +568,7 @@ def vapor_mass_transfer_eqn(blk, t, x, j): self.vapor_phase.length_domain, initialize=100, units=lunits("power") / lunits("temperature") / lunits("length"), - doc="Vapor-liquid heat transfer coefficient", + doc="Vapor-liquid heat transfer coefficient multiplied by heat transfer area per unit column length", ) # Heat transfer @@ -664,19 +674,34 @@ def pressure_at_interface(blk, t, x, j): else: zb = self.liquid_phase.length_domain.prev(x) lprops = blk.liquid_phase.properties[t, zb] - return blk.pressure_equil[t, x, j] == (lprops.fug_phase_comp["Liq", j]) + return blk.mass_transfer_driving_force[t, x, j] == ( + blk.vapor_phase.properties[t, x].mole_frac_comp[j] + * pyunits.convert( + blk.vapor_phase.properties[t, x].pressure, + to_units=lunits("pressure"), + ) + - lprops.fug_phase_comp["Liq", j] + ) # ========================================================================= # Scaling routine def calculate_scaling_factors(self): super().calculate_scaling_factors() + vunits = ( + self.config.vapor_phase.property_package.get_metadata().get_derived_units + ) + lunits = ( + self.config.liquid_phase.property_package.get_metadata().get_derived_units + ) + # --------------------------------------------------------------------- # Scale variables - for (t, x, j), v in self.pressure_equil.items(): + # TODO revisit + for (t, x, j), v in self.mass_transfer_driving_force.items(): if iscale.get_scaling_factor(v) is None: sf_pe = iscale.get_scaling_factor( - self.pressure_equil, default=None, warning=True + self.mass_transfer_driving_force, default=None, warning=True ) if sf_pe is None: sf_pe = iscale.get_scaling_factor( @@ -691,7 +716,7 @@ def calculate_scaling_factors(self): warning=True, ) - iscale.set_scaling_factor(v, sf_pe) + iscale.set_scaling_factor(v, sf_pe*20) for (t, x), v in self.vapor_phase.heat.items(): if iscale.get_scaling_factor(v) is None: @@ -721,7 +746,7 @@ def calculate_scaling_factors(self): iscale.constraint_scaling_transform( v, iscale.get_scaling_factor( - self.pressure_equil[t, x, j], default=1, warning=False + self.mass_transfer_driving_force[t, x, j], default=1, warning=False ), ) @@ -734,13 +759,16 @@ def calculate_scaling_factors(self): ) for (t, x, j), v in self.liquid_mass_transfer_eqn.items(): + zf = self.vapor_phase.length_domain.next(x) try: sf = iscale.get_scaling_factor( - self.interphase_mass_transfer[t, x, j], default=1, warning=False + self.interphase_mass_transfer[t, zf, j], default=1, warning=False ) except KeyError: # This implies a non-volatile component - sf = 1 + sf = iscale.get_scaling_factor( + self.liquid_phase.mass_transfer_term[t, x, "Liq", j], default=1, warning=True + ) iscale.constraint_scaling_transform(v, sf) for (t, x, j), v in self.vapor_mass_transfer_eqn.items(): @@ -748,9 +776,18 @@ def calculate_scaling_factors(self): sf = iscale.get_scaling_factor( self.interphase_mass_transfer[t, x, j], default=1, warning=False ) + # Account for the fact that this equation is written on a vapor unit basis + sf_units = pyunits.convert_value( + 1, + from_units=1 / (lunits("amount") / lunits("time") / lunits("length")), + to_units=1 / (vunits("amount") / vunits("time") / vunits("length")) + ) + sf *= sf_units except KeyError: # This implies a non-volatile component - sf = 1 + sf = iscale.get_scaling_factor( + self.vapor_phase.mass_transfer_term[t, x, "Vap", j], default=1, warning=True + ) iscale.constraint_scaling_transform(v, sf) for (t, x), v in self.heat_transfer_eqn1.items(): diff --git a/idaes/models_extra/column_models/solvent_condenser.py b/idaes/models_extra/column_models/solvent_condenser.py index 465c080793..972da62d09 100644 --- a/idaes/models_extra/column_models/solvent_condenser.py +++ b/idaes/models_extra/column_models/solvent_condenser.py @@ -61,7 +61,7 @@ class SolventCondenserData(UnitModelBlockData): """ Condenser unit for solvent column models using separate property packages - for liquid and vpor phases. + for liquid and vapor phases. Unit model to condense the vapor from the top of a solvent column. """ @@ -321,11 +321,11 @@ def build(self): flow_basis = self.liquid_phase[t_init].get_material_flow_basis() if flow_basis == MaterialFlowBasis.molar: fb = "flow_mole" - elif flow_basis == MaterialFlowBasis.molar: - fb = "flow_mass" + # elif flow_basis == MaterialFlowBasis.mass: + # fb = "flow_mass" else: raise ConfigurationError( - f"{self.name} SolventCondenser only supports mass or molar " + f"{self.name} SolventCondenser only supports molar " f"basis for MaterialFlowBasis." ) @@ -450,7 +450,14 @@ def calculate_scaling_factors(self): ), ) elif j in self.liquid_phase.component_list: - iscale.constraint_scaling_transform(v, value(1 / self.zero_flow_param)) + iscale.constraint_scaling_transform( + v, + iscale.get_scaling_factor( + self.liquid_phase[t].get_material_flow_terms("Liq",j), + default=1, + warning=True, + ), + ) else: pass # no need to scale this constraint