diff --git a/setup.py b/setup.py index 750b967f..d8565034 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ python_requires=">=3.8", install_requires=[ # "watertap @ https://github.com/watertap-org/watertap/archive/main.zip", # uncomment if we need to point to main mid release cycle - "watertap>=1.0.0rc0", + "watertap==1.0.0rc0", "idaes-pse==2.5.0", "pyomo==6.7.3", "nrel-pysam == 5.1.0", diff --git a/src/watertap_contrib/reflo/costing/units/multi_effect_crystallizer.py b/src/watertap_contrib/reflo/costing/units/multi_effect_crystallizer.py new file mode 100644 index 00000000..f5191ebe --- /dev/null +++ b/src/watertap_contrib/reflo/costing/units/multi_effect_crystallizer.py @@ -0,0 +1,358 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + +import pyomo.environ as pyo +from idaes.core.util.misc import StrEnum +from idaes.core.util.constants import Constants +from idaes.core.util.exceptions import ConfigurationError +from watertap.costing.util import register_costing_parameter_block +from watertap_contrib.reflo.costing.util import make_capital_cost_var + + +class MultiEffectCrystallizerCostType(StrEnum): + mass_basis = "mass_basis" + volume_basis = "volume_basis" + + +def build_recovered_nacl_cost_param_block(blk): + + blk.cost = pyo.Param( + mutable=True, + initialize=0, + doc="Recovered cost (sale price) of NaCl", + units=pyo.units.USD_2020 / pyo.units.kg, + ) + + costing = blk.parent_block() + costing.register_flow_type("NaCl_recovered", blk.cost) + + +def build_steam_cost_param_block(blk): + + blk.cost = pyo.Param( + mutable=True, + initialize=0.004, + doc="Steam cost", + units=pyo.units.USD_2018 * pyo.units.meter**-3, + ) + + costing = blk.parent_block() + costing.register_flow_type("steam", blk.cost) + + +def build_multi_effect_crystallizer_cost_param_block(blk): + + # Crystallizer operating cost information from literature + blk.fob_unit_cost = pyo.Var( + initialize=675000, + doc="Forced circulation crystallizer reference free-on-board cost (Woods, 2007)", + units=pyo.units.USD_2007, + ) + + blk.ref_capacity = pyo.Var( + initialize=1, + doc="Forced circulation crystallizer reference crystal capacity (Woods, 2007)", + units=pyo.units.kg / pyo.units.s, + ) + + blk.ref_exponent = pyo.Var( + initialize=0.53, + doc="Forced circulation crystallizer cost exponent factor (Woods, 2007)", + units=pyo.units.dimensionless, + ) + + blk.iec_percent = pyo.Var( + initialize=1.43, + doc="Forced circulation crystallizer installed equipment cost (Diab and Gerogiorgis, 2017)", + units=pyo.units.dimensionless, + ) + + blk.volume_cost = pyo.Var( + initialize=16320, + doc="Forced circulation crystallizer cost per volume (Yusuf et al., 2019)", + units=pyo.units.USD_2007, ## TODO: Needs confirmation, but data is from Perry apparently + ) + + blk.vol_basis_exponent = pyo.Var( + initialize=0.47, + doc="Forced circulation crystallizer volume-based cost exponent (Yusuf et al., 2019)", + units=pyo.units.dimensionless, + ) + + blk.heat_exchanger_capital_factor = pyo.Var( + initialize=420, + units=pyo.units.USD_2018 / (pyo.units.meter**2), + doc="Heat exchanger cost per area", # TODO: need reference and costing year + ) + + blk.heat_exchanger_endplates_capital_factor = pyo.Var( + initialize=1020, + units=pyo.units.USD_2018, + doc="Heat exchanger endplates cost per area", # TODO: need reference and costing year + ) + + blk.heat_exchanger_endplates_capital_basis = pyo.Var( + initialize=10, + units=pyo.units.meter**2, + doc="Heat exchanger endplates cost per area basis", # TODO: need reference and costing year + ) + + blk.heat_exchanger_endplates_capital_exponent = pyo.Var( + initialize=0.6, + units=pyo.units.dimensionless, + doc="Heat exchanger endplates cost per area basis", # TODO: need reference and costing year + ) + + +@register_costing_parameter_block( + build_rule=build_recovered_nacl_cost_param_block, + parameter_block_name="nacl_recovered", +) +@register_costing_parameter_block( + build_rule=build_steam_cost_param_block, + parameter_block_name="steam", +) +@register_costing_parameter_block( + build_rule=build_multi_effect_crystallizer_cost_param_block, + parameter_block_name="multi_effect_crystallizer", +) +def cost_multi_effect_crystallizer( + blk, cost_type=MultiEffectCrystallizerCostType.mass_basis +): + """ + Function for costing the forced circulation crystallizer by the mass flow of produced crystals. + The operating cost model assumes that heat is supplied via condensation of saturated steam (see Dutta et al.) + + Args: + cost_type: Option for crystallizer cost function type - volume or mass basis + """ + global costing_package + + costing_package = blk.costing_package + make_capital_cost_var(blk) + costing_package.add_cost_factor(blk, "TIC") + + if cost_type == MultiEffectCrystallizerCostType.mass_basis: + effect_costing_method = cost_crystallizer_effect_by_crystal_mass + elif cost_type == MultiEffectCrystallizerCostType.volume_basis: + effect_costing_method = cost_crystallizer_effect_by_volume + else: + raise ConfigurationError( + f"{blk.unit_model.name} received invalid argument for cost_type:" + f" {cost_type}. Argument must be a member of the MultiEffectCrystallizerCostType Enum." + ) + + total_capex_expr = 0 + + if not hasattr(blk.unit_model, "effects"): + + assert blk.unit_model.config.standalone + # blk.unit_model is CrystallizerEffect as standalone unit + effect_capex_expr = 0 + + # Add capital of crystallizer + effect_capex_expr += effect_costing_method(blk.unit_model) + # separate expression for crystallizer capex for reporting + blk.add_component( + f"capital_cost_crystallizer", + pyo.Expression( + expr=blk.cost_factor * effect_costing_method(blk.unit_model), + doc=f"Capital cost for only crystallizer including TIC", + ), + ) + + # Add capital of heat exchangers + effect_capex_expr += cost_crystallizer_heat_exchanger(blk.unit_model) + # separate expression for heat exchanger for reporting + blk.add_component( + f"capital_cost_heat_exchanger", + pyo.Expression( + expr=blk.cost_factor * cost_crystallizer_heat_exchanger(blk.unit_model), + doc=f"Capital cost for only heat exchanger including TIC", + ), + ) + + blk.capital_cost_constraint = pyo.Constraint( + expr=blk.capital_cost == blk.cost_factor * effect_capex_expr, + doc=f"Constraint for total capital cost of crystallizer.", + ) + + _cost_effect_flows(blk.unit_model, 1) + + elif hasattr(blk.unit_model, "effects"): + + # blk.unit_model is MultiEffectCrystallizer + + for effect_number, eff in blk.unit_model.effects.items(): + + effect_capex_expr = 0 + + effect_capex_var = pyo.Var( + initialize=1e5, + units=costing_package.base_currency, + doc=f"Capital cost effect {effect_number}", + ) + + # Add capital of crystallizer + effect_capex_expr += effect_costing_method(eff.effect) + # separate expression for crystallizer capex for reporting + blk.add_component( + f"capital_cost_crystallizer_effect_{effect_number}", + pyo.Expression( + expr=blk.cost_factor * effect_costing_method(eff.effect), + doc=f"Capital cost for only crystallizer for effect {effect_number} including TIC", + ), + ) + + # Add capital of heat exchangers + effect_capex_expr += cost_crystallizer_heat_exchanger(eff.effect) + # separate expression for heat exchanger for reporting + blk.add_component( + f"capital_cost_heat_exchanger_effect_{effect_number}", + pyo.Expression( + expr=blk.cost_factor * cost_crystallizer_heat_exchanger(eff.effect), + doc=f"Capital cost for only heat exchanger for effect {effect_number} including TIC", + ), + ) + + # Total effect CAPEX constraint + effect_capex_constr = pyo.Constraint( + expr=effect_capex_var == effect_capex_expr, + doc=f"Constraint for total capital cost of effect {effect_number}.", + ) + # Add effect Var and Constraint to costing blk + blk.add_component(f"capital_cost_effect_{effect_number}", effect_capex_var) + blk.add_component( + f"total_capital_cost_effect_{effect_number}_constraint", + effect_capex_constr, + ) + total_capex_expr += effect_capex_expr + _cost_effect_flows(eff.effect, effect_number) + + blk.capital_cost_constraint = pyo.Constraint( + expr=blk.capital_cost == blk.cost_factor * total_capex_expr + ) + else: + raise ConfigurationError( + f"{blk.unit_model.name} is not a valid unit model type for this costing package:" + f" {type(blk.unit_model)}. Must be either CrystallizerEffect or MultiEffectCrystallizer." + ) + + +def cost_crystallizer_heat_exchanger(effect): + + capital_cost_hx_effect = 0 + + capital_cost_hx = pyo.units.convert( + costing_package.multi_effect_crystallizer.heat_exchanger_capital_factor + * effect.heat_exchanger_area, + to_units=costing_package.base_currency, + ) + + capital_cost_hx_effect += capital_cost_hx + + dimensionless_hx_area = pyo.units.convert( + ( + effect.heat_exchanger_area + / costing_package.multi_effect_crystallizer.heat_exchanger_endplates_capital_basis + ), + to_units=pyo.units.dimensionless, + ) + + capital_cost_hx_endplates = pyo.units.convert( + costing_package.multi_effect_crystallizer.heat_exchanger_endplates_capital_factor + * (dimensionless_hx_area) + ** costing_package.multi_effect_crystallizer.heat_exchanger_endplates_capital_exponent, + to_units=costing_package.base_currency, + ) + + capital_cost_hx_effect += capital_cost_hx_endplates + + return capital_cost_hx_effect + + +def cost_crystallizer_effect_by_crystal_mass(effect): + """ + Mass-based capital cost for forced circulation crystallizer + """ + + capital_cost_effect = pyo.units.convert( + ( + costing_package.multi_effect_crystallizer.iec_percent + * costing_package.multi_effect_crystallizer.fob_unit_cost + * ( + sum( + effect.properties_solids[0].flow_mass_phase_comp["Sol", j] + for j in effect.config.property_package.solute_set + ) + / costing_package.multi_effect_crystallizer.ref_capacity + ) + ** costing_package.multi_effect_crystallizer.ref_exponent + ), + to_units=costing_package.base_currency, + ) + + return capital_cost_effect + + +def cost_crystallizer_effect_by_volume(effect): + """ + Volume-based capital cost for forced circulation crystallizer + """ + + capital_cost_effect = pyo.units.convert( + ( + costing_package.multi_effect_crystallizer.volume_cost + * ( + ( + pyo.units.convert( + effect.volume_suspension + * (effect.height_crystallizer / effect.height_slurry), + to_units=(pyo.units.ft) ** 3, + ) + ) + / pyo.units.ft**3 + ) + ** costing_package.multi_effect_crystallizer.vol_basis_exponent + ), + to_units=costing_package.base_currency, + ) + + return capital_cost_effect + + +def _cost_effect_flows(effect, effect_number): + + costing_package.cost_flow( + pyo.units.convert( + ( + effect.magma_circulation_flow_vol + * effect.dens_mass_slurry + * Constants.acceleration_gravity + * effect.pump_head_height + / effect.efficiency_pump + ), + to_units=pyo.units.kW, + ), + "electricity", + ) + + costing_package.cost_flow( + effect.properties_solids[0].flow_mass_phase_comp["Sol", "NaCl"], + "NaCl_recovered", + ) + + if effect_number == 1: + costing_package.cost_flow( + effect.heating_steam[0].flow_vol_phase["Vap"], + "steam", + ) diff --git a/src/watertap_contrib/reflo/unit_models/crystallizer_effect.py b/src/watertap_contrib/reflo/unit_models/crystallizer_effect.py new file mode 100644 index 00000000..da3956ca --- /dev/null +++ b/src/watertap_contrib/reflo/unit_models/crystallizer_effect.py @@ -0,0 +1,429 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + +from copy import deepcopy + +# Import Pyomo libraries +from pyomo.environ import ( + Var, + Param, + check_optimal_termination, + units as pyunits, +) +from pyomo.common.config import ConfigValue + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + useDefault, +) +from idaes.core.util.config import is_physical_parameter_block +from idaes.core.util.exceptions import InitializationError +import idaes.core.util.scaling as iscale +import idaes.logger as idaeslog + +from watertap.core.solvers import get_solver +from watertap.core.util.initialization import interval_initializer +from watertap.unit_models.crystallizer import CrystallizationData +from watertap.unit_models.mvc.components.lmtd_chen_callback import ( + delta_temperature_chen_callback, +) + +from watertap_contrib.reflo.costing.units.multi_effect_crystallizer import ( + cost_multi_effect_crystallizer, +) + +__author__ = "Oluwamayowa Amusat, Zhuoran Zhang, Kurban Sitterley" + + +@declare_process_block_class("CrystallizerEffect") +class CrystallizerEffectData(CrystallizationData): + """ + Zero dimensional model for crystallizer effect + """ + + CONFIG = CrystallizationData.CONFIG() + + CONFIG.declare( + "property_package_vapor", + ConfigValue( + default=useDefault, + domain=is_physical_parameter_block, + description="Property package to use for heating steam properties", + doc="""Property parameter object used to define steam property calculations, + **default** - useDefault. + **Valid values:** { + **useDefault** - use default package from parent model or flowsheet, + **PhysicalParameterObject** - a PhysicalParameterBlock object.}""", + ), + ) + CONFIG.declare( + "standalone", + ConfigValue( + default=True, + domain=bool, + description="Flag to indicate if model is used alone or as part of MultiEffectCrystallizer.", + doc="""Flag to indicate if model is used alone or as part of MultiEffectCrystallizer unit model, + **default** - True.""", + ), + ) + + def build(self): + super().build() + + tmp_dict = dict(**self.config.property_package_args) + tmp_dict["has_phase_equilibrium"] = False + tmp_dict["parameters"] = self.config.property_package + tmp_dict["defined_state"] = False + + self.properties_pure_water = self.config.property_package.state_block_class( + self.flowsheet().config.time, + doc="Material properties of pure water vapour at outlet", + **tmp_dict, + ) + + self.add_port(name="pure_water", block=self.properties_pure_water) + + self.properties_pure_water[0].flow_mass_phase_comp["Vap", "H2O"].fix(0) + self.properties_pure_water[0].flow_mass_phase_comp["Liq", "NaCl"].fix(0) + self.properties_pure_water[0].flow_mass_phase_comp["Sol", "NaCl"].fix(0) + self.properties_in[0].conc_mass_phase_comp["Liq", "NaCl"] + self.properties_in[0].flow_vol_phase["Liq"] + + self.inlet.temperature.setub(1000) + self.outlet.temperature.setub(1000) + self.solids.temperature.setub(1000) + self.vapor.temperature.setub(1000) + self.pure_water.temperature.setub(1000) + + self.steam_pressure = Param( + initialize=3, + mutable=True, + units=pyunits.bar, + doc="Steam pressure (gauge) for crystallizer heating: 3 bar default based on Dutta example", + ) + + self.efficiency_pump = Param( + initialize=0.7, + mutable=True, + units=pyunits.dimensionless, + doc="Crystallizer pump efficiency - assumed", + ) + + self.pump_head_height = Param( + initialize=1, + mutable=True, + units=pyunits.m, + doc="Crystallizer pump head height - assumed, unvalidated", + ) + + self.energy_flow_superheated_vapor = Var( + initialize=1e5, + bounds=(-5e6, 5e6), + units=pyunits.kilowatt, + doc="Energy that could be supplied from vapor", + ) + + self.delta_temperature_in = Var( + self.flowsheet().time, + initialize=35, + bounds=(None, None), + units=pyunits.degK, + doc="Temperature difference at the inlet side", + ) + + self.delta_temperature_out = Var( + self.flowsheet().time, + initialize=35, + bounds=(None, None), + units=pyunits.degK, + doc="Temperature difference at the outlet side", + ) + + delta_temperature_chen_callback(self) + + self.heat_exchanger_area = Var( + initialize=1000.0, + bounds=(0, None), + units=pyunits.m**2, + doc="Heat exchanger area", + ) + + self.overall_heat_transfer_coefficient = Var( + initialize=0.1, + bounds=(0, None), + units=pyunits.kilowatt / pyunits.m**2 / pyunits.degK, + doc="Overall heat transfer coefficient for heat exchangers", + ) + + @self.Constraint() + def eq_pure_water_mass_flow_rate(b): + return ( + b.properties_pure_water[0].flow_mass_phase_comp["Liq", "H2O"] + == b.properties_vapor[0].flow_mass_phase_comp["Vap", "H2O"] + ) + + @self.Constraint(doc="Thermal energy in the vapor") + def eq_vapor_energy_constraint(b): + return b.energy_flow_superheated_vapor == ( + b.properties_vapor[0].flow_mass_phase_comp["Vap", "H2O"] + * ( + b.properties_pure_water[ + 0 + ].dh_vap_mass_solvent # Latent heat from the vapor + + b.properties_vapor[0].enth_mass_solvent["Vap"] + - b.properties_pure_water[0].enth_mass_solvent["Vap"] + ) + ) + + self.del_component(self.eq_p_con1) + self.del_component(self.eq_p_con2) + + @self.Constraint() + def eq_p_con1(b): + return b.pressure_operating == b.properties_out[0].pressure + + @self.Constraint() + def eq_p_con2(b): + return b.pressure_operating == b.properties_solids[0].pressure + + @self.Constraint() + def eq_p_con4(b): + return b.properties_pure_water[0].pressure == self.pressure_operating + + @self.Constraint() + def eq_p_con5(b): + return b.properties_pure_water[0].pressure_sat == self.pressure_operating + + if self.config.standalone: + + tmp_dict["parameters"] = self.config.property_package_vapor + tmp_dict["defined_state"] = False + + self.heating_steam = self.config.property_package_vapor.state_block_class( + self.flowsheet().config.time, + doc="Material properties of inlet heating steam", + **tmp_dict, + ) + + self.add_port(name="steam", block=self.heating_steam) + self.steam.temperature.setub(1000) + + @self.Constraint(doc="Change in temperature at inlet") + def eq_delta_temperature_inlet(b): + return ( + b.delta_temperature_in[0] + == b.heating_steam[0].temperature - b.temperature_operating + ) + + @self.Constraint(doc="Change in temperature at outlet") + def eq_delta_temperature_outlet(b): + return ( + b.delta_temperature_out[0] + == b.heating_steam[0].temperature - b.properties_in[0].temperature + ) + + @self.Constraint(doc="Heat transfer equation") + def eq_heat_transfer(b): + return b.work_mechanical[0] == ( + b.overall_heat_transfer_coefficient + * b.heat_exchanger_area + * b.delta_temperature[0] + ) + + @self.Constraint(doc="Calculate mass flow rate of heating steam") + def eq_heating_steam_flow_rate(b): + return b.work_mechanical[0] == ( + pyunits.convert( + b.heating_steam[0].dh_vap_mass + * b.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"], + to_units=pyunits.kJ / pyunits.s, + ) + ) + + else: + self.del_component(self.inlet) + self.del_component(self.outlet) + self.del_component(self.solids) + self.del_component(self.vapor) + self.del_component(self.pure_water) + + def initialize_build( + self, + state_args=None, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + ): + """ + General wrapper for pressure changer initialization routines + + Keyword Arguments: + state_args : a dict of arguments to be passed to the property + package(s) to provide an initial state for + initialization (see documentation of the specific + property package) (default = {}). + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None) + solver : str indicating which solver to use during + initialization (default = None) + + Returns: None + """ + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") + + opt = get_solver(solver, optarg) + + flags = self.properties_in.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args, + hold_state=True, + ) + init_log.info_high("Initialization Step 1 Complete.") + + if state_args is None: + state_args = {} + state_dict = self.properties_in[ + self.flowsheet().config.time.first() + ].define_port_members() + + for k in state_dict.keys(): + if state_dict[k].is_indexed(): + state_args[k] = {} + for m in state_dict[k].keys(): + state_args[k][m] = state_dict[k][m].value + else: + state_args[k] = state_dict[k].value + + self.properties_out.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args, + ) + + state_args_solids = deepcopy(state_args) + + for p, j in self.properties_solids.phase_component_set: + if p == "Sol": + state_args_solids["flow_mass_phase_comp"][p, j] = state_args[ + "flow_mass_phase_comp" + ]["Liq", j] + elif p in ["Liq", "Vap"]: + state_args_solids["flow_mass_phase_comp"][p, j] = 1e-8 + + self.properties_solids.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args_solids, + ) + + state_args_vapor = deepcopy(state_args) + + for p, j in self.properties_vapor.phase_component_set: + if p == "Vap": + state_args_vapor["flow_mass_phase_comp"][p, j] = state_args[ + "flow_mass_phase_comp" + ]["Liq", j] + elif p == "Liq" or p == "Sol": + state_args_vapor["flow_mass_phase_comp"][p, j] = 1e-8 + + self.properties_vapor.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args_vapor, + ) + + self.properties_pure_water.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args_vapor, + ) + + if hasattr(self, "heating_steam"): + state_args_steam = deepcopy(state_args) + + for p, j in self.properties_vapor.phase_component_set: + state_args_steam["flow_mass_phase_comp"][p, j] = 1 + + self.heating_steam.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args_steam, + ) + + init_log.info_high("Initialization Step 2 Complete.") + + interval_initializer(self) + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res))) + + self.properties_in.release_state(flags, outlvl=outlvl) + init_log.info_high( + "Initialization Complete: {}".format(idaeslog.condition(res)) + ) + + if not check_optimal_termination(res): + raise InitializationError(f"Unit model {self.name} failed to initialize") + + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + if iscale.get_scaling_factor(self.work_mechanical[0]) is None: + iscale.set_scaling_factor( + self.work_mechanical[0], + iscale.get_scaling_factor( + self.properties_in[0].flow_mass_phase_comp["Vap", "H2O"] + ) + * iscale.get_scaling_factor( + self.properties_in[0].enth_mass_solvent["Vap"] + ) + * 1e3, + ) + if iscale.get_scaling_factor(self.energy_flow_superheated_vapor) is None: + iscale.set_scaling_factor( + self.energy_flow_superheated_vapor, + iscale.get_scaling_factor( + self.properties_vapor[0].flow_mass_phase_comp["Vap", "H2O"] + ) + * iscale.get_scaling_factor( + self.properties_vapor[0].enth_mass_solvent["Vap"] + ), + ) + + if iscale.get_scaling_factor(self.delta_temperature_in[0]) is None: + iscale.set_scaling_factor(self.delta_temperature_in[0], 0.1) + + if iscale.get_scaling_factor(self.delta_temperature_out[0]) is None: + iscale.set_scaling_factor(self.delta_temperature_out[0], 0.1) + + if iscale.get_scaling_factor(self.heat_exchanger_area) is None: + iscale.set_scaling_factor(self.heat_exchanger_area, 0.1) + + if iscale.get_scaling_factor(self.overall_heat_transfer_coefficient) is None: + iscale.set_scaling_factor(self.overall_heat_transfer_coefficient, 0.01) + + for _, c in self.eq_p_con4.items(): + sf = iscale.get_scaling_factor(self.properties_pure_water[0].pressure) + iscale.constraint_scaling_transform(c, sf) + + @property + def default_costing_method(self): + return cost_multi_effect_crystallizer diff --git a/src/watertap_contrib/reflo/unit_models/multi_effect_crystallizer.py b/src/watertap_contrib/reflo/unit_models/multi_effect_crystallizer.py new file mode 100644 index 00000000..97e60b06 --- /dev/null +++ b/src/watertap_contrib/reflo/unit_models/multi_effect_crystallizer.py @@ -0,0 +1,491 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + + +# Import Pyomo libraries +from pyomo.environ import ( + ConcreteModel, + check_optimal_termination, + assert_optimal_termination, + Var, + Constraint, + Expression, + Suffix, + RangeSet, + units as pyunits, +) +from pyomo.common.config import ConfigBlock, ConfigValue, In, PositiveInt + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + MaterialBalanceType, + UnitModelBlockData, + useDefault, + FlowsheetBlock, + UnitModelCostingBlock, +) +from idaes.core.util.exceptions import InitializationError, ConfigurationError +import idaes.core.util.scaling as iscale +from idaes.core.util.config import is_physical_parameter_block +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.tables import create_stream_table_dataframe +import idaes.logger as idaeslog +from idaes.core.util.constants import Constants + +from watertap.core import InitializationMixin, ControlVolume0DBlock +from watertap.core.solvers import get_solver + +from watertap_contrib.reflo.unit_models.crystallizer_effect import CrystallizerEffect +from watertap_contrib.reflo.costing.units.multi_effect_crystallizer import ( + cost_multi_effect_crystallizer, +) + +_log = idaeslog.getLogger(__name__) + +__author__ = "Oluwamayowa Amusat, Zhuoran Zhang, Kurban Sitterley" + + +@declare_process_block_class("MultiEffectCrystallizer") +class MultiEffectCrystallizerData(InitializationMixin, UnitModelBlockData): + + 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. The filtration unit does not support dynamic + behavior, thus this must be False.""", + ), + ) + + CONFIG.declare( + "has_holdup", + ConfigValue( + default=False, + domain=In([False]), + description="Holdup construction flag - must be False", + doc="""Indicates whether holdup terms should be constructed or not. + **default** - False. The filtration unit does not have defined volume, thus + this must be False.""", + ), + ) + + CONFIG.declare( + "property_package", + ConfigValue( + default=useDefault, + domain=is_physical_parameter_block, + description="Property package to use for control volume", + doc="""Property parameter object used to define property calculations, + **default** - useDefault. + **Valid values:** { + **useDefault** - use default package from parent model or flowsheet, + **PhysicalParameterObject** - a PhysicalParameterBlock object.}""", + ), + ) + + CONFIG.declare( + "property_package_vapor", + ConfigValue( + default=useDefault, + domain=is_physical_parameter_block, + description="Property package to use for heating steam properties", + doc="""Property parameter object used to define steam property calculations, + **default** - useDefault. + **Valid values:** { + **useDefault** - use default package from parent model or flowsheet, + **PhysicalParameterObject** - a PhysicalParameterBlock 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( + "number_effects", + ConfigValue( + default=4, + domain=PositiveInt, + description="Number of effects of the multi-effect crystallizer system", + doc="""A ConfigBlock specifying the number of effects, which can only be 4.""", + ), + ) + + def build(self): + super().build() + + if self.config.number_effects <= 1: + raise ConfigurationError( + "The MultiEffectCrystallizer model requires more than 1 effect." + "To model a crystallizer with one effect, use the CrystallizerEffect model with 'standalone=True'." + ) + + self.scaling_factor = Suffix(direction=Suffix.EXPORT) + + self.control_volume = ControlVolume0DBlock( + dynamic=False, + has_holdup=False, + property_package=self.config.property_package, + property_package_args=self.config.property_package_args, + ) + + self.control_volume.add_state_blocks(has_phase_equilibrium=False) + self.control_volume.add_material_balances( + balance_type=MaterialBalanceType.componentPhase, has_mass_transfer=True + ) + self.add_inlet_port(name="inlet", block=self.control_volume) + self.add_outlet_port(name="outlet", block=self.control_volume) + + self.number_effects = self.config.number_effects + self.Effects = RangeSet(self.config.number_effects) + + self.first_effect = self.Effects.first() + self.last_effect = self.Effects.last() + + self.effects = FlowsheetBlock(self.Effects, dynamic=False) + + # There is no solid NaCl coming in + self.control_volume.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(0) + # All solid NaCl accessed through properties_solids + self.control_volume.mass_transfer_term[0, "Sol", "NaCl"].fix(0) + + # Local expressions to aggregate material flows for all effects + total_mass_flow_water_in_expr = 0 + total_mass_flow_salt_in_expr = 0 + total_mass_flow_pure_water_out_expr = 0 + total_mass_flow_water_out_expr = 0 + total_mass_flow_salt_out_expr = 0 + total_flow_vol_in_expr = 0 + total_flow_vol_out_expr = 0 + + for n, eff in self.effects.items(): + eff.effect = effect = CrystallizerEffect( + property_package=self.config.property_package, + property_package_vapor=self.config.property_package_vapor, + standalone=False, + ) + effect.properties_in[0].conc_mass_phase_comp + + total_flow_vol_in_expr += effect.properties_in[0].flow_vol_phase["Liq"] + total_flow_vol_out_expr += effect.properties_pure_water[0].flow_vol_phase[ + "Liq" + ] + + total_mass_flow_water_in_expr += effect.properties_in[ + 0 + ].flow_mass_phase_comp["Liq", "H2O"] + + total_mass_flow_salt_in_expr += effect.properties_in[ + 0 + ].flow_mass_phase_comp["Liq", "NaCl"] + + total_mass_flow_pure_water_out_expr += effect.properties_pure_water[ + 0 + ].flow_mass_phase_comp["Liq", "H2O"] + + total_mass_flow_water_out_expr += effect.properties_out[ + 0 + ].flow_mass_phase_comp["Liq", "H2O"] + + total_mass_flow_salt_out_expr += ( + effect.properties_out[0].flow_mass_phase_comp["Liq", "NaCl"] + + effect.properties_solids[0].flow_mass_phase_comp["Sol", "NaCl"] + ) + + if n == self.first_effect: + tmp_dict = dict(**self.config.property_package_args) + tmp_dict["has_phase_equilibrium"] = False + tmp_dict["parameters"] = self.config.property_package_vapor + tmp_dict["defined_state"] = False + + effect.heating_steam = ( + self.config.property_package_vapor.state_block_class( + self.flowsheet().config.time, + doc="Material properties of inlet heating steam", + **tmp_dict, + ) + ) + + @effect.Constraint( + doc="Change in temperature at inlet for first effect" + ) + def eq_delta_temperature_inlet_effect_1(b): + return ( + b.delta_temperature_in[0] + == b.heating_steam[0].temperature - b.temperature_operating + ) + + @effect.Constraint( + doc="Change in temperature at outlet for first effect" + ) + def eq_delta_temperature_outlet_effect_1(b): + return ( + b.delta_temperature_out[0] + == b.heating_steam[0].temperature + - b.properties_in[0].temperature + ) + + @effect.Constraint(doc="Heat transfer equation for first effect") + def eq_heat_transfer_effect_1(b): + return ( + b.work_mechanical[0] + == b.overall_heat_transfer_coefficient + * b.heat_exchanger_area + * b.delta_temperature[0] + ) + + @effect.Constraint(doc="Calculate mass flow rate of heating steam") + def eq_heating_steam_flow_rate(b): + return b.work_mechanical[0] == ( + pyunits.convert( + b.heating_steam[0].dh_vap_mass + * b.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"], + to_units=pyunits.kJ / pyunits.s, + ) + ) + + self.add_port(name="solids", block=effect.properties_solids) + self.add_port(name="vapor", block=effect.properties_vapor) + self.add_port(name="steam", block=effect.heating_steam) + self.steam.temperature.setub(1000) + + else: + + prev_effect = self.effects[n - 1].effect + + del_temp_in_constr = Constraint( + expr=effect.delta_temperature_in[0] + == prev_effect.properties_vapor[0].temperature + - effect.temperature_operating, + doc=f"Change in temperature at inlet for effect {n}", + ) + effect.add_component( + f"eq_delta_temperature_inlet_effect_{n}", del_temp_in_constr + ) + + del_temp_out_constr = Constraint( + expr=effect.delta_temperature_out[0] + == prev_effect.properties_pure_water[0].temperature + - effect.properties_in[0].temperature, + doc=f"Change in temperature at outlet for effect {n}", + ) + effect.add_component( + f"eq_delta_temperature_outlet_effect_{n}", del_temp_out_constr + ) + + hx_constr = Constraint( + expr=prev_effect.energy_flow_superheated_vapor + == effect.overall_heat_transfer_coefficient + * effect.heat_exchanger_area + * effect.delta_temperature[0], + doc=f"Heat transfer equation for effect {n}", + ) + effect.add_component(f"eq_heat_transfer_effect_{n}", hx_constr) + + energy_flow_constr = Constraint( + expr=effect.work_mechanical[0] + == prev_effect.energy_flow_superheated_vapor, + doc=f"Energy supplied to effect {n}", + ) + effect.add_component( + f"eq_energy_for_effect_{n}_from_effect_{n - 1}", energy_flow_constr + ) + + @self.Constraint(doc="Mass transfer term for liquid water") + def eq_mass_transfer_term_liq_water(b): + return b.control_volume.mass_transfer_term[0, "Liq", "H2O"] == -1 * ( + total_mass_flow_water_out_expr + ) + + @self.Constraint(doc="Mass transfer term for vapor water") + def eq_mass_transfer_term_vap_water(b): + return b.control_volume.mass_transfer_term[0, "Vap", "H2O"] == -1 * ( + b.effects[1].effect.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"] + ) + + @self.Constraint(doc="Mass transfer term for salt in liquid phase") + def eq_mass_transfer_term_liq_salt(b): + return b.control_volume.mass_transfer_term[0, "Liq", "NaCl"] == -1 * ( + total_mass_flow_salt_out_expr + ) + + @self.Constraint(doc="Steam flow") + def eq_overall_steam_flow(b): + return ( + b.control_volume.properties_in[0].flow_mass_phase_comp["Vap", "H2O"] + == b.effects[1] + .effect.heating_steam[0] + .flow_mass_phase_comp["Vap", "H2O"] + ) + + @self.Constraint(doc="Mass balance of water for all effects") + def eq_overall_mass_balance_water_in(b): + return ( + b.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == total_mass_flow_water_in_expr + ) + + @self.Constraint(doc="Mass balance of salt for all effects") + def eq_overall_mass_balance_salt_in(b): + return ( + b.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + == total_mass_flow_salt_in_expr + ) + + @self.Constraint(doc="Control volume temperature at outlet") + def eq_temperature_outlet(b): + return ( + b.control_volume.properties_out[0].temperature + == b.effects[b.last_effect].effect.properties_pure_water[0].temperature + ) + + @self.Constraint(doc="Control volume pressure at outlet") + def eq_isobaric(b): + return ( + b.control_volume.properties_out[0].pressure + == b.control_volume.properties_in[0].pressure + ) + + self.total_flow_vol_in = Expression(expr=total_flow_vol_in_expr) + + self.recovery_vol_phase = Var( + ["Liq"], + initialize=0.75, + bounds=(0, 1), + units=pyunits.dimensionless, + doc="Overall unit recovery on volumetric basis", + ) + + @self.Constraint(doc="Volumetric recovery") + def eq_recovery_vol_phase(b): + return ( + b.recovery_vol_phase["Liq"] + == total_flow_vol_out_expr / total_flow_vol_in_expr + ) + + def initialize_build( + self, + state_args=None, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + ): + """ + Keyword Arguments: + state_args : a dict of arguments to be passed to the property + package(s) to provide an initial state for + initialization (see documentation of the specific + property package) (default = {}). + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None) + solver : str indicating which solver to use during + initialization (default = None) + + Returns: None + """ + + init_args = dict( + state_args=state_args, + outlvl=outlvl, + solver=solver, + optarg=optarg, + ) + + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") + + flow_mass_water_initial = ( + self.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value + / self.number_effects + ) + flow_mass_salt_initial = ( + self.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .value + / self.number_effects + ) + + opt = get_solver(solver, optarg) + + for n, eff in self.effects.items(): + # Each effect is first solved in a vacuum with linking constraints deactivated + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].fix( + flow_mass_water_initial + ) + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].fix( + flow_mass_salt_initial + ) + if n == 1: + if not degrees_of_freedom(eff.effect) == 0: + raise InitializationError( + f"Degrees of freedom in first effect must be zero during initialization, " + f"but has {degrees_of_freedom(eff.effect)}. Check inlet conditions and re-initialize." + ) + eff.effect.initialize(**init_args) + inlet_conc = ( + eff.effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + .value + ) + mass_transfer_coeff = eff.effect.overall_heat_transfer_coefficient.value + else: + # Deactivate contraint that links energy flow between effects + linking_constr = getattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + linking_constr.deactivate() + eff.effect.initialize(**init_args) + linking_constr.activate() + # Fix inlet feed concentration to be equal across all effects + eff.effect.properties_in[0].conc_mass_phase_comp["Liq", "NaCl"].fix( + inlet_conc + ) + eff.effect.overall_heat_transfer_coefficient.fix(mass_transfer_coeff) + + # Unfix inlet mass flow rates to allow unit model to determine based on energy flows + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].unfix() + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].unfix() + + init_log.info(f"Initialization of Effect {n} Complete.") + + with idaeslog.solver_log(init_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + init_log.info("Initialization Step 3 {}.".format(idaeslog.condition(res))) + + if not check_optimal_termination(res): + raise InitializationError(f"Unit model {self.name} failed to initialize") + + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + + if iscale.get_scaling_factor(self.recovery_vol_phase) is None: + iscale.set_scaling_factor(self.recovery_vol_phase, 10) + + @property + def default_costing_method(self): + return cost_multi_effect_crystallizer diff --git a/src/watertap_contrib/reflo/unit_models/tests/test_crystallizer_effect.py b/src/watertap_contrib/reflo/unit_models/tests/test_crystallizer_effect.py new file mode 100644 index 00000000..83baade7 --- /dev/null +++ b/src/watertap_contrib/reflo/unit_models/tests/test_crystallizer_effect.py @@ -0,0 +1,447 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + +import pytest + +from pyomo.environ import ( + ConcreteModel, + Var, + Param, + value, + assert_optimal_termination, + units as pyunits, +) +from pyomo.network import Port + +from idaes.core import FlowsheetBlock, UnitModelCostingBlock, MaterialFlowBasis +from idaes.core.util.testing import initialization_tester +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_variables, + number_total_constraints, + number_unused_variables, +) +from idaes.core.util.scaling import ( + calculate_scaling_factors, + unscaled_variables_generator, + badly_scaled_var_generator, +) +from pyomo.util.check_units import assert_units_consistent +import idaes.logger as idaeslog + +from watertap.core.solvers import get_solver +from watertap.property_models.unit_specific.cryst_prop_pack import ( + NaClParameterBlock, + NaClParameterData, + NaClStateBlock, +) +from watertap.property_models.water_prop_pack import ( + WaterParameterBlock, + WaterStateBlock, +) + +from watertap_contrib.reflo.unit_models.crystallizer_effect import CrystallizerEffect +from watertap_contrib.reflo.costing import TreatmentCosting + +# Get default solver for testing +solver = get_solver() + + +def build_effect(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = NaClParameterBlock() + m.fs.vapor_properties = WaterParameterBlock() + + m.fs.unit = eff = CrystallizerEffect( + property_package=m.fs.properties, + property_package_vapor=m.fs.vapor_properties, + standalone=True, + ) + + feed_flow_mass = 1 + feed_mass_frac_NaCl = 0.15 + atm_pressure = 101325 * pyunits.Pa + saturated_steam_pressure_gage = 3 * pyunits.bar + saturated_steam_pressure = atm_pressure + pyunits.convert( + saturated_steam_pressure_gage, to_units=pyunits.Pa + ) + feed_temperature = 273.15 + 20 + crystallizer_yield = 0.5 + operating_pressure_eff = 0.45 + + feed_mass_frac_H2O = 1 - feed_mass_frac_NaCl + eps = 1e-6 + + eff.inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].fix( + feed_flow_mass * feed_mass_frac_NaCl + ) + eff.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + eff.inlet.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(eps) + eff.inlet.flow_mass_phase_comp[0, "Vap", "H2O"].fix(eps) + + eff.inlet.pressure[0].fix(atm_pressure) + eff.inlet.temperature[0].fix(feed_temperature) + + eff.heating_steam[0].pressure_sat + + eff.heating_steam.calculate_state( + var_args={ + ("flow_mass_phase_comp", ("Liq", "H2O")): 0, + ("pressure", None): saturated_steam_pressure, + ("pressure_sat", None): saturated_steam_pressure, + }, + hold_state=True, + ) + eff.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"].unfix() + + eff.crystallization_yield["NaCl"].fix(crystallizer_yield) + eff.crystal_growth_rate.fix() + eff.souders_brown_constant.fix() + eff.crystal_median_length.fix() + eff.overall_heat_transfer_coefficient.fix(0.1) + + eff.pressure_operating.fix(operating_pressure_eff * pyunits.bar) + + return m + + +class TestCrystallizerEffect: + @pytest.fixture(scope="class") + def effect_frame(self): + m = build_effect() + return m + + @pytest.mark.unit + def test_config(self, effect_frame): + m = effect_frame + + assert len(m.fs.unit.config) == 6 + + assert not m.fs.unit.config.dynamic + assert not m.fs.unit.config.has_holdup + assert m.fs.unit.config.property_package is m.fs.properties + assert m.fs.unit.config.property_package_vapor is m.fs.vapor_properties + assert m.fs.unit.config.standalone + assert_units_consistent(m) + + @pytest.mark.unit + def test_build(self, effect_frame): + m = effect_frame + + # test ports and variables + port_lst = ["inlet", "outlet", "solids", "vapor", "steam", "pure_water"] + port_vars_lst = ["flow_mass_phase_comp", "pressure", "temperature"] + state_blks = [ + "properties_in", + "properties_out", + "properties_pure_water", + "properties_solids", + "properties_vapor", + "heating_steam", + ] + + for port_str in port_lst: + assert hasattr(m.fs.unit, port_str) + port = getattr(m.fs.unit, port_str) + assert len(port.vars) == 3 + assert isinstance(port, Port) + for var_str in port_vars_lst: + assert hasattr(port, var_str) + var = getattr(port, var_str) + assert isinstance(var, Var) + + for b in state_blks: + assert hasattr(m.fs.unit, b) + sb = getattr(m.fs.unit, b) + assert sb[0].temperature.ub == 1000 + if b == "heating_steam": + assert isinstance(sb, WaterStateBlock) + else: + assert isinstance(sb, NaClStateBlock) + + effect_vars = [ + "crystal_median_length", + "crystal_growth_rate", + "souders_brown_constant", + "crystallization_yield", + "product_volumetric_solids_fraction", + "temperature_operating", + "pressure_operating", + "dens_mass_magma", + "dens_mass_slurry", + "work_mechanical", + "diameter_crystallizer", + "height_slurry", + "height_crystallizer", + "magma_circulation_flow_vol", + "relative_supersaturation", + "energy_flow_superheated_vapor", + "delta_temperature_in", + "delta_temperature_out", + "heat_exchanger_area", + "overall_heat_transfer_coefficient", + ] + + for ev in effect_vars: + v = getattr(m.fs.unit, ev) + assert isinstance(v, Var) + + effect_params = [ + "approach_temperature_heat_exchanger", + "dimensionless_crystal_length", + "steam_pressure", + "efficiency_pump", + "pump_head_height", + ] + + for ep in effect_params: + p = getattr(m.fs.unit, ep) + assert isinstance(p, Param) + + assert number_variables(m) == 288 + assert number_total_constraints(m) == 120 + assert number_unused_variables(m) == 44 + + assert_units_consistent(m) + + @pytest.mark.unit + def test_dof(self, effect_frame): + m = effect_frame + assert degrees_of_freedom(m) == 0 + + @pytest.mark.unit + def test_calculate_scaling(self, effect_frame): + m = effect_frame + + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Liq", "NaCl") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Vap", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Sol", "NaCl") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Vap", "H2O") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + + calculate_scaling_factors(m) + unscaled_var_list = list(unscaled_variables_generator(m)) + assert len(unscaled_var_list) == 0 + badly_scaled_var_list = list(badly_scaled_var_generator(m)) + assert len(badly_scaled_var_list) == 0 + + @pytest.mark.component + def test_initialize(self, effect_frame): + m = effect_frame + initialization_tester(m) + + @pytest.mark.component + def test_solve(self, effect_frame): + m = effect_frame + results = solver.solve(m) + assert_optimal_termination(results) + + @pytest.mark.unit + def test_conservation(self, effect_frame): + + m = effect_frame + + comp_lst = ["NaCl", "H2O"] + phase_lst = ["Sol", "Liq", "Vap"] + + phase_comp_list = [ + (p, j) + for j in comp_lst + for p in phase_lst + if (p, j) in m.fs.unit.properties_in[0].phase_component_set + ] + flow_mass_in = sum( + m.fs.unit.properties_in[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_out = sum( + m.fs.unit.properties_out[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_solids = sum( + m.fs.unit.properties_solids[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_vapor = sum( + m.fs.unit.properties_vapor[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + + assert ( + abs( + value(flow_mass_in - flow_mass_out - flow_mass_solids - flow_mass_vapor) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + flow_mass_in * m.fs.unit.properties_in[0].enth_mass_phase["Liq"] + - flow_mass_out * m.fs.unit.properties_out[0].enth_mass_phase["Liq"] + - flow_mass_vapor + * m.fs.unit.properties_vapor[0].enth_mass_solvent["Vap"] + - flow_mass_solids + * m.fs.unit.properties_solids[0].enth_mass_solute["Sol"] + - flow_mass_solids + * m.fs.unit.properties_solids[0].dh_crystallization_mass_comp[ + "NaCl" + ] + + m.fs.unit.work_mechanical[0] + ) + ) + <= 1e-2 + ) + + @pytest.mark.unit + def test_solution(self, effect_frame): + m = effect_frame + + eff_dict = { + "product_volumetric_solids_fraction": 0.132975, + "temperature_operating": 359.48, + "dens_mass_magma": 281.24, + "dens_mass_slurry": 1298.23, + "work_mechanical": {0.0: 1704.47}, + "diameter_crystallizer": 1.0799, + "height_slurry": 1.073, + "height_crystallizer": 1.883, + "magma_circulation_flow_vol": 0.119395, + "relative_supersaturation": {"NaCl": 0.567038}, + "t_res": 1.0228, + "volume_suspension": 0.982965, + "eq_max_allowable_velocity": 2.6304, + "eq_vapor_space_height": 0.809971, + "eq_minimum_height_diameter_ratio": 1.6199, + "energy_flow_superheated_vapor": 1518.73, + "delta_temperature_in": {0.0: 57.39}, + "delta_temperature_out": {0.0: 123.72}, + "delta_temperature": {0.0: 86.31}, + "heat_exchanger_area": 197.47, + } + + for v, r in eff_dict.items(): + effv = getattr(m.fs.unit, v) + if effv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(effv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(effv), rel=1e-3) == r + + steam_dict = { + "flow_mass_phase_comp": {("Liq", "H2O"): 0.0, ("Vap", "H2O"): 0.799053}, + "temperature": 416.87, + "pressure": 401325.0, + "dh_vap_mass": 2133119.0, + "pressure_sat": 401325.0, + } + + for v, r in steam_dict.items(): + sv = getattr(m.fs.unit.heating_steam[0], v) + if sv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(sv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(sv), rel=1e-3) == r + + @pytest.mark.component + def test_costing(self, effect_frame): + m = effect_frame + m.fs.costing = TreatmentCosting() + m.fs.unit.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + ) + m.fs.costing.nacl_recovered.cost.set_value(-0.024) + m.fs.costing.cost_process() + m.fs.costing.add_LCOW(m.fs.unit.properties_in[0].flow_vol_phase["Liq"]) + m.fs.costing.add_specific_energy_consumption( + m.fs.unit.properties_in[0].flow_vol_phase["Liq"], name="SEC" + ) + results = solver.solve(m) + assert_optimal_termination(results) + + sys_costing_dict = { + "aggregate_capital_cost": 868242.67, + "aggregate_flow_electricity": 2.1715, + "aggregate_flow_NaCl_recovered": 0.075, + "aggregate_flow_steam": 0.383078, + "aggregate_flow_costs": { + "electricity": 1564.26, + "NaCl_recovered": -67456.42, + "steam": 56766.9, + }, + "aggregate_direct_capital_cost": 434121.33, + "total_capital_cost": 868242.67, + "total_operating_cost": 16922.02, + "maintenance_labor_chemical_operating_cost": 26047.28, + "total_fixed_operating_cost": 26047.28, + "total_variable_operating_cost": -9125.25, + "total_annualized_cost": 114126.95, + "LCOW": 4.0094, + "SEC": 0.66875, + } + + for v, r in sys_costing_dict.items(): + cv = getattr(m.fs.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + eff_costing_dict = { + "capital_cost": 868242.67, + "direct_capital_cost": 434121.33, + "capital_cost_crystallizer": 659171.48, + "capital_cost_heat_exchanger": 209071.19, + } + + for v, r in eff_costing_dict.items(): + cv = getattr(m.fs.unit.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + assert ( + pytest.approx( + value(m.fs.unit.heating_steam[0].flow_vol_phase["Vap"]), rel=1e-3 + ) + == 0.383078 + ) + + assert value(m.fs.unit.heating_steam[0].flow_vol_phase["Liq"]) == 0 diff --git a/src/watertap_contrib/reflo/unit_models/tests/test_multi_effect_crystallizer.py b/src/watertap_contrib/reflo/unit_models/tests/test_multi_effect_crystallizer.py new file mode 100644 index 00000000..8044c076 --- /dev/null +++ b/src/watertap_contrib/reflo/unit_models/tests/test_multi_effect_crystallizer.py @@ -0,0 +1,2288 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + +import pytest + +from pyomo.util.check_units import assert_units_consistent +from pyomo.environ import ( + ConcreteModel, + Var, + Objective, + assert_optimal_termination, + value, + units as pyunits, +) +from pyomo.network import Port + +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_variables, + number_total_constraints, + number_unused_variables, +) +from idaes.core.util.scaling import ( + calculate_scaling_factors, + unscaled_variables_generator, + badly_scaled_var_generator, +) +from idaes.core import UnitModelCostingBlock +from idaes.core.util.exceptions import ConfigurationError + +from watertap.core.solvers import get_solver +from watertap.property_models.unit_specific.cryst_prop_pack import ( + NaClParameterBlock, + NaClStateBlock, +) +from watertap.property_models.water_prop_pack import ( + WaterParameterBlock, + WaterStateBlock, +) + +from watertap_contrib.reflo.costing import TreatmentCosting +from watertap_contrib.reflo.unit_models.multi_effect_crystallizer import ( + MultiEffectCrystallizer, +) +from watertap_contrib.reflo.unit_models.crystallizer_effect import CrystallizerEffect + +__author__ = "Kurban Sitterley" + +solver = get_solver() +feed_pressure = 101325 +feed_temperature = 273.15 + 20 +eps = 1e-12 + + +def build_mec2(): + + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = NaClParameterBlock() + m.fs.vapor_properties = WaterParameterBlock() + + m.fs.unit = mec = MultiEffectCrystallizer( + property_package=m.fs.properties, + property_package_vapor=m.fs.vapor_properties, + number_effects=2, + ) + + num_effects = 2 + feed_flow_mass = 2 + total_feed_flow_mass = num_effects * feed_flow_mass + feed_mass_frac_NaCl = 0.45 + crystallizer_yield = 0.6 + operating_pressures = [0.45, 0.25] + feed_mass_frac_H2O = 1 - feed_mass_frac_NaCl + + atm_pressure = 101325 * pyunits.Pa + saturated_steam_pressure_gage = 5 * pyunits.bar + saturated_steam_pressure = atm_pressure + pyunits.convert( + saturated_steam_pressure_gage, to_units=pyunits.Pa + ) + + for (_, eff), op_pressure in zip(mec.effects.items(), operating_pressures): + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].set_value( + feed_flow_mass * feed_mass_frac_H2O + ) + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].set_value( + feed_flow_mass * feed_mass_frac_NaCl + ) + eff.effect.properties_in[0].pressure.fix(feed_pressure) + eff.effect.properties_in[0].temperature.fix(feed_temperature) + + eff.effect.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(eps) + eff.effect.properties_in[0].flow_mass_phase_comp["Vap", "H2O"].fix(eps) + eff.effect.properties_in[0].conc_mass_phase_comp[...] + eff.effect.crystallization_yield["NaCl"].fix(crystallizer_yield) + eff.effect.crystal_growth_rate.fix() + eff.effect.souders_brown_constant.fix() + eff.effect.crystal_median_length.fix() + eff.effect.pressure_operating.fix( + pyunits.convert(op_pressure * pyunits.bar, to_units=pyunits.Pa) + ) + eff.effect.overall_heat_transfer_coefficient.set_value(0.1) + + first_effect = m.fs.unit.effects[1].effect + first_effect.overall_heat_transfer_coefficient.fix(0.1) + first_effect.heating_steam[0].dh_vap_mass + first_effect.heating_steam.calculate_state( + var_args={ + ("flow_mass_phase_comp", ("Liq", "H2O")): 0, + ("pressure", None): saturated_steam_pressure, + ("pressure_sat", None): saturated_steam_pressure, + }, + hold_state=True, + ) + first_effect.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"].unfix() + + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].fix( + total_feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].fix( + total_feed_flow_mass * feed_mass_frac_NaCl + ) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(0) + m.fs.unit.control_volume.properties_in[0].pressure.fix(feed_pressure) + m.fs.unit.control_volume.properties_in[0].temperature.fix(feed_temperature) + + return m + + +def build_mec3(): + + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = NaClParameterBlock() + m.fs.vapor_properties = WaterParameterBlock() + + m.fs.unit = mec = MultiEffectCrystallizer( + property_package=m.fs.properties, + property_package_vapor=m.fs.vapor_properties, + number_effects=3, + ) + + num_effects = 3 + feed_flow_mass = 10 + total_feed_flow_mass = num_effects * feed_flow_mass + feed_mass_frac_NaCl = 0.25 + crystallizer_yield = 0.55 + operating_pressures = [0.45, 0.25, 0.208] + feed_mass_frac_H2O = 1 - feed_mass_frac_NaCl + + atm_pressure = 101325 * pyunits.Pa + saturated_steam_pressure_gage = 3.5 * pyunits.bar + saturated_steam_pressure = atm_pressure + pyunits.convert( + saturated_steam_pressure_gage, to_units=pyunits.Pa + ) + + for (_, eff), op_pressure in zip(mec.effects.items(), operating_pressures): + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].fix( + feed_flow_mass * feed_mass_frac_NaCl + ) + eff.effect.properties_in[0].pressure.fix(feed_pressure) + eff.effect.properties_in[0].temperature.fix(feed_temperature) + + eff.effect.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(eps) + eff.effect.properties_in[0].flow_mass_phase_comp["Vap", "H2O"].fix(eps) + eff.effect.properties_in[0].conc_mass_phase_comp[...] + eff.effect.crystallization_yield["NaCl"].fix(crystallizer_yield) + eff.effect.crystal_growth_rate.fix() + eff.effect.souders_brown_constant.fix() + eff.effect.crystal_median_length.fix() + eff.effect.pressure_operating.fix( + pyunits.convert(op_pressure * pyunits.bar, to_units=pyunits.Pa) + ) + eff.effect.overall_heat_transfer_coefficient.set_value(0.1) + + first_effect = m.fs.unit.effects[1].effect + first_effect.overall_heat_transfer_coefficient.fix(0.1) + first_effect.heating_steam[0].dh_vap_mass + first_effect.heating_steam.calculate_state( + var_args={ + ("flow_mass_phase_comp", ("Liq", "H2O")): 0, + ("pressure", None): saturated_steam_pressure, + ("pressure_sat", None): saturated_steam_pressure, + }, + hold_state=True, + ) + first_effect.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"].unfix() + + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].fix( + total_feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].fix( + total_feed_flow_mass * feed_mass_frac_NaCl + ) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(0) + m.fs.unit.control_volume.properties_in[0].pressure.fix(feed_pressure) + m.fs.unit.control_volume.properties_in[0].temperature.fix(feed_temperature) + + return m + + +def build_mec4(): + + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = NaClParameterBlock() + m.fs.vapor_properties = WaterParameterBlock() + + m.fs.unit = mec = MultiEffectCrystallizer( + property_package=m.fs.properties, property_package_vapor=m.fs.vapor_properties + ) + + operating_pressures = [0.45, 0.25, 0.208, 0.095] + + feed_flow_mass = 1 + feed_mass_frac_NaCl = 0.15 + crystallizer_yield = 0.5 + feed_mass_frac_H2O = 1 - feed_mass_frac_NaCl + + total_feed_flow_mass = 4 + + atm_pressure = 101325 * pyunits.Pa + saturated_steam_pressure_gage = 3 * pyunits.bar + saturated_steam_pressure = atm_pressure + pyunits.convert( + saturated_steam_pressure_gage, to_units=pyunits.Pa + ) + + for (_, eff), op_pressure in zip(mec.effects.items(), operating_pressures): + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].fix( + feed_flow_mass * feed_mass_frac_H2O + ) + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].fix( + feed_flow_mass * feed_mass_frac_NaCl + ) + eff.effect.properties_in[0].pressure.fix(feed_pressure) + eff.effect.properties_in[0].temperature.fix(feed_temperature) + + eff.effect.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(eps) + eff.effect.properties_in[0].flow_mass_phase_comp["Vap", "H2O"].fix(eps) + eff.effect.properties_in[0].conc_mass_phase_comp[...] + eff.effect.crystallization_yield["NaCl"].fix(crystallizer_yield) + eff.effect.crystal_growth_rate.fix() + eff.effect.souders_brown_constant.fix() + eff.effect.crystal_median_length.fix() + eff.effect.pressure_operating.fix( + pyunits.convert(op_pressure * pyunits.bar, to_units=pyunits.Pa) + ) + eff.effect.overall_heat_transfer_coefficient.set_value(0.1) + + first_effect = m.fs.unit.effects[1].effect + + first_effect.overall_heat_transfer_coefficient.fix(0.1) + first_effect.heating_steam[0].pressure_sat + first_effect.heating_steam[0].dh_vap_mass + first_effect.heating_steam.calculate_state( + var_args={ + ("flow_mass_phase_comp", ("Liq", "H2O")): 0, + ("pressure", None): saturated_steam_pressure, + ("pressure_sat", None): saturated_steam_pressure, + }, + hold_state=True, + ) + first_effect.heating_steam[0].flow_mass_phase_comp["Vap", "H2O"].unfix() + + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].fix( + total_feed_flow_mass * feed_mass_frac_H2O + ) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].fix( + total_feed_flow_mass * feed_mass_frac_NaCl + ) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp["Sol", "NaCl"].fix(0) + m.fs.unit.control_volume.properties_in[0].pressure.fix(feed_pressure) + m.fs.unit.control_volume.properties_in[0].temperature.fix(feed_temperature) + + return m + + +@pytest.mark.component +def test_single_effect_error(): + + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.properties = NaClParameterBlock() + m.fs.vapor_properties = WaterParameterBlock() + + error_msg = ( + "The MultiEffectCrystallizer model requires more than 1 effect." + "To model a crystallizer with one effect, use the CrystallizerEffect model with 'standalone=True'." + ) + + with pytest.raises(ConfigurationError, match=error_msg): + m.fs.unit = MultiEffectCrystallizer( + property_package=m.fs.properties, + property_package_vapor=m.fs.vapor_properties, + number_effects=1, + ) + + +class TestMultiEffectCrystallizer_2Effects: + @pytest.fixture(scope="class") + def MEC2_frame(self): + """ + Test crystallizer with 2 effects + """ + m = build_mec2() + return m + + @pytest.mark.unit + def test_config(self, MEC2_frame): + m = MEC2_frame + assert len(m.fs.unit.config) == 6 + + assert not m.fs.unit.config.dynamic + assert not m.fs.unit.config.has_holdup + assert m.fs.unit.config.property_package is m.fs.properties + assert m.fs.unit.config.property_package_vapor is m.fs.vapor_properties + assert m.fs.unit.config.number_effects == 2 + assert m.fs.unit.config.number_effects == len(m.fs.unit.effects) + + for _, eff in m.fs.unit.effects.items(): + assert isinstance(eff.effect, CrystallizerEffect) + assert not eff.effect.config.standalone + + @pytest.mark.unit + def test_build(self, MEC2_frame): + m = MEC2_frame + + # test ports and variables + port_lst = ["inlet", "outlet", "solids", "vapor", "steam"] + port_vars_lst = ["flow_mass_phase_comp", "pressure", "temperature"] + state_blks = [ + "properties_in", + "properties_out", + "properties_pure_water", + "properties_solids", + "properties_vapor", + ] + + for port_str in port_lst: + assert hasattr(m.fs.unit, port_str) + port = getattr(m.fs.unit, port_str) + assert len(port.vars) == 3 + assert isinstance(port, Port) + for var_str in port_vars_lst: + assert hasattr(port, var_str) + var = getattr(port, var_str) + assert isinstance(var, Var) + + for n, eff in m.fs.unit.effects.items(): + for b in state_blks: + assert hasattr(eff.effect, b) + sb = getattr(eff.effect, b) + assert isinstance(sb, NaClStateBlock) + assert sb[0].temperature.ub == 1000 + + effect_params = [ + "approach_temperature_heat_exchanger", + "dimensionless_crystal_length", + ] + effect_vars = [ + "crystal_growth_rate", + "crystal_median_length", + "crystallization_yield", + "dens_mass_magma", + "dens_mass_slurry", + "diameter_crystallizer", + "height_crystallizer", + "height_slurry", + "magma_circulation_flow_vol", + "pressure_operating", + "product_volumetric_solids_fraction", + "relative_supersaturation", + "souders_brown_constant", + "t_res", + "temperature_operating", + "volume_suspension", + "work_mechanical", + ] + + effect_exprs = [ + "delta_temperature", + "eq_max_allowable_velocity", + "eq_minimum_height_diameter_ratio", + "eq_vapor_space_height", + ] + + effect_constr = [ + "eq_mass_balance_constraints", + "eq_solubility_massfrac_equality_constraint", + "eq_removal_balance", + "eq_vol_fraction_solids", + "eq_dens_magma", + "eq_operating_pressure_constraint", + "eq_relative_supersaturation", + "eq_enthalpy_balance", + "eq_p_con3", + "eq_T_con1", + "eq_T_con2", + "eq_T_con3", + "eq_minimum_hex_circulation_rate_constraint", + "eq_dens_mass_slurry", + "eq_residence_time", + "eq_suspension_volume", + "eq_vapor_head_diameter_constraint", + "eq_slurry_height_constraint", + "eq_crystallizer_height_constraint", + "eq_pure_water_mass_flow_rate", + "eq_vapor_energy_constraint", + "eq_p_con1", + "eq_p_con2", + "eq_p_con4", + "eq_p_con5", + ] + + assert hasattr(m.fs.unit, "recovery_vol_phase") + assert isinstance(m.fs.unit.recovery_vol_phase, Var) + + for n, eff in m.fs.unit.effects.items(): + for p in effect_params: + assert hasattr(eff.effect, p) + for v in effect_vars: + assert hasattr(eff.effect, v) + for e in effect_exprs: + assert hasattr(eff.effect, e) + for c in effect_constr: + assert hasattr(eff.effect, c) + assert hasattr(eff.effect, f"eq_delta_temperature_inlet_effect_{n}") + assert hasattr(eff.effect, f"eq_delta_temperature_outlet_effect_{n}") + assert hasattr(eff.effect, f"eq_heat_transfer_effect_{n}") + if n == 1: + assert number_variables(eff.effect) == 154 + assert number_total_constraints(eff.effect) == 128 + assert number_unused_variables(eff.effect) == 2 + assert hasattr(eff.effect, "heating_steam") + assert isinstance(eff.effect.heating_steam, WaterStateBlock) + assert eff.effect.heating_steam[0].temperature.ub == 1000 + assert hasattr(eff.effect, "eq_heating_steam_flow_rate") + if n != 1: + assert number_variables(eff.effect) == 148 + assert number_total_constraints(eff.effect) == 126 + assert number_unused_variables(eff.effect) == 0 + assert hasattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + + assert number_variables(m) == 461 + assert number_total_constraints(m) == 267 + assert number_unused_variables(m) == 43 + + assert_units_consistent(m) + + @pytest.mark.unit + def test_dof(self, MEC2_frame): + m = MEC2_frame + # With mass flow rates into each of the effects unfixed, + # the model is under specified + assert degrees_of_freedom(m) == 2 + + # Fixing flow rates into individual effects will reduce DOF... + for n, eff in m.fs.unit.effects.items(): + eff.effect.properties_in[0].flow_mass_phase_comp.fix() + if n == 1: + assert degrees_of_freedom(eff.effect) == 0 + else: + assert degrees_of_freedom(eff.effect) == 3 + # ... and result in an overspecified model. + assert degrees_of_freedom(m) == -2 + + for n, eff in m.fs.unit.effects.items(): + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].unfix() + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].unfix() + assert degrees_of_freedom(m) == 2 + + @pytest.mark.unit + def test_calculate_scaling(self, MEC2_frame): + m = MEC2_frame + + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "NaCl") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 10, index=("Vap", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Sol", "NaCl") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Vap", "H2O") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + + calculate_scaling_factors(m) + unscaled_var_list = list(unscaled_variables_generator(m)) + assert len(unscaled_var_list) == 0 + badly_scaled_var_list = list(badly_scaled_var_generator(m)) + assert len(badly_scaled_var_list) == 0 + + @pytest.mark.component + def test_initialize(self, MEC2_frame): + m = MEC2_frame + m.fs.unit.initialize() + c0 = value( + m.fs.unit.effects[1] + .effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + ) + htc = value(m.fs.unit.effects[1].effect.overall_heat_transfer_coefficient) + assert degrees_of_freedom(m) == 0 + for n, eff in m.fs.unit.effects.items(): + assert ( + not eff.effect.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .is_fixed() + ) + assert ( + not eff.effect.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .is_fixed() + ) + assert ( + pytest.approx( + value( + eff.effect.properties_in[0].conc_mass_phase_comp["Liq", "NaCl"] + ), + rel=1e-6, + ) + == c0 + ) + if n == 1: + assert degrees_of_freedom(eff.effect) == 2 + if n != 1: + assert ( + eff.effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + .is_fixed() + ) + linking_constr = getattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + assert linking_constr.active + assert eff.effect.overall_heat_transfer_coefficient.is_fixed() + assert value(eff.effect.overall_heat_transfer_coefficient) == htc + assert degrees_of_freedom(eff.effect) == 3 + + @pytest.mark.component + def test_solve(self, MEC2_frame): + + m = MEC2_frame + + results = solver.solve(m) + assert_optimal_termination(results) + + @pytest.mark.component + def test_conservation(self, MEC2_frame): + + m = MEC2_frame + + comp_lst = ["NaCl", "H2O"] + phase_lst = ["Sol", "Liq", "Vap"] + + total_mass_flow_water_in = 0 + total_mass_flow_salt_in = 0 + total_mass_flow_water_out = 0 + + for _, e in m.fs.unit.effects.items(): + eff = e.effect + + phase_comp_list = [ + (p, j) + for j in comp_lst + for p in phase_lst + if (p, j) in eff.properties_in[0].phase_component_set + ] + flow_mass_in = sum( + eff.properties_in[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_out = sum( + eff.properties_out[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_solids = sum( + eff.properties_solids[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_vapor = sum( + eff.properties_vapor[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + + assert ( + abs( + value( + flow_mass_in + - flow_mass_out + - flow_mass_solids + - flow_mass_vapor + ) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + flow_mass_in * eff.properties_in[0].enth_mass_phase["Liq"] + - flow_mass_out * eff.properties_out[0].enth_mass_phase["Liq"] + - flow_mass_vapor + * eff.properties_vapor[0].enth_mass_solvent["Vap"] + - flow_mass_solids + * eff.properties_solids[0].enth_mass_solute["Sol"] + - flow_mass_solids + * eff.properties_solids[0].dh_crystallization_mass_comp["NaCl"] + + eff.work_mechanical[0] + ) + ) + <= 1e-2 + ) + + total_mass_flow_water_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + ) + total_mass_flow_salt_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + total_mass_flow_water_out += value( + eff.properties_pure_water[0].flow_mass_phase_comp["Liq", "H2O"] + ) + + # Test control volume mass balance + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .value, + rel=1e-6, + ) + == total_mass_flow_salt_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_out[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_out + ) + + @pytest.mark.component + def test_solution(self, MEC2_frame): + m = MEC2_frame + + assert ( + pytest.approx(value(m.fs.unit.recovery_vol_phase["Liq"]), rel=1e-3) + == 0.10289 + ) + + unit_results_dict = { + 1: { + "temperature_operating": 359.4, + "pressure_operating": 45000.0, + "dens_mass_magma": 395.317, + "dens_mass_slurry": 1349.047, + "work_mechanical": {0.0: 370.3}, + "diameter_crystallizer": 0.50292, + "height_slurry": 22.853, + "height_crystallizer": 23.23, + "magma_circulation_flow_vol": 0.026733, + "relative_supersaturation": {"NaCl": 0.76747}, + "t_res": 1.0228, + "volume_suspension": 4.5398, + "eq_max_allowable_velocity": 2.630, + "eq_vapor_space_height": 0.37719, + "eq_minimum_height_diameter_ratio": 0.75438, + "energy_flow_superheated_vapor": 329.356, + "delta_temperature_in": {0.0: 72.57}, + "delta_temperature_out": {0.0: 138.9}, + "delta_temperature": {0.0: 102.1}, + "heat_exchanger_area": 36.253, + }, + 2: { + "temperature_operating": 344.8, + "pressure_operating": 25000.0, + "dens_mass_magma": 392.7, + "dens_mass_slurry": 1352.0, + "work_mechanical": {0.0: 329.3}, + "diameter_crystallizer": 0.6011, + "height_slurry": 19.582, + "height_crystallizer": 20.03, + "magma_circulation_flow_vol": 0.02368, + "relative_supersaturation": {"NaCl": 0.7739}, + "t_res": 1.0228, + "volume_suspension": 5.558, + "eq_max_allowable_velocity": 3.464, + "eq_vapor_space_height": 0.45086, + "eq_minimum_height_diameter_ratio": 0.9017, + "energy_flow_superheated_vapor": 364.13, + "delta_temperature_in": {0.0: 14.61}, + "delta_temperature_out": {0.0: 58.77}, + "delta_temperature": {0.0: 31.59}, + "heat_exchanger_area": 104.251, + }, + } + + for n, d in unit_results_dict.items(): + eff = m.fs.unit.effects[n].effect + for v, r in d.items(): + effv = getattr(eff, v) + if effv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(effv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(effv), rel=1e-3) == r + + steam_results_dict = { + "flow_mass_phase_comp": {("Vap", "H2O"): 0.177583}, + "temperature": 432.06, + "pressure": 601325.0, + "dh_vap_mass": 2085530.7, + "pressure_sat": 601325.0, + } + + for v, r in steam_results_dict.items(): + sv = getattr(m.fs.unit.effects[1].effect.heating_steam[0], v) + if sv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(sv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(sv), rel=1e-3) == r + + @pytest.mark.component + def test_costing(self, MEC2_frame): + m = MEC2_frame + m.fs.costing = TreatmentCosting() + # m.fs.costing.base_currency = pyunits.USD_2018 + m.fs.unit.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method_arguments={"cost_type": "mass_basis"}, + ) + + m.fs.costing.nacl_recovered.cost.set_value(-0.011) + m.fs.costing.cost_process() + m.fs.costing.add_LCOW(m.fs.unit.total_flow_vol_in) + m.fs.costing.add_specific_energy_consumption( + m.fs.unit.total_flow_vol_in, name="SEC" + ) + results = solver.solve(m) + assert_optimal_termination(results) + + sys_costing_dict = { + "aggregate_capital_cost": 3902338.7, + "aggregate_flow_electricity": 0.95390, + "aggregate_flow_NaCl_recovered": 1.0799, + "aggregate_flow_steam": 0.05888, + "aggregate_flow_costs": { + "electricity": 687.144, + "NaCl_recovered": -445206.6, + "steam": 8726.6, + }, + "aggregate_direct_capital_cost": 1951169.3, + "total_capital_cost": 3902338.7, + "total_operating_cost": -318722.72, + "maintenance_labor_chemical_operating_cost": 117070.16, + "total_fixed_operating_cost": 117070.16, + "total_variable_operating_cost": -435792.88, + "total_annualized_cost": 118167.30, + "LCOW": 1.27058636, + "SEC": 0.08991097, + } + + for v, r in sys_costing_dict.items(): + cv = getattr(m.fs.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + eff_costing_dict = { + "capital_cost": 3902338.7, + "direct_capital_cost": 1951169.3, + "capital_cost_crystallizer_effect_1": 1777273.3, + "capital_cost_heat_exchanger_effect_1": 40936.8, + "capital_cost_effect_1": 909105.0, + "capital_cost_crystallizer_effect_2": 1971550.7, + "capital_cost_heat_exchanger_effect_2": 112577.7, + "capital_cost_effect_2": 1042064.2, + } + + for v, r in eff_costing_dict.items(): + cv = getattr(m.fs.unit.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + +class TestMultiEffectCrystallizer_3Effects: + @pytest.fixture(scope="class") + def MEC3_frame(self): + """ + Test 3 effect crystallizer + """ + m = build_mec3() + return m + + @pytest.mark.unit + def test_config(self, MEC3_frame): + m = MEC3_frame + assert len(m.fs.unit.config) == 6 + + assert not m.fs.unit.config.dynamic + assert not m.fs.unit.config.has_holdup + assert m.fs.unit.config.property_package is m.fs.properties + assert m.fs.unit.config.property_package_vapor is m.fs.vapor_properties + assert m.fs.unit.config.number_effects == 3 + assert m.fs.unit.config.number_effects == len(m.fs.unit.effects) + + assert isinstance(m.fs.unit.effects, FlowsheetBlock) + + for _, eff in m.fs.unit.effects.items(): + assert isinstance(eff.effect, CrystallizerEffect) + assert not eff.effect.config.standalone + + @pytest.mark.unit + def test_build(self, MEC3_frame): + m = MEC3_frame + + # test ports and variables + port_lst = ["inlet", "outlet", "solids", "vapor", "steam"] + port_vars_lst = ["flow_mass_phase_comp", "pressure", "temperature"] + state_blks = [ + "properties_in", + "properties_out", + "properties_pure_water", + "properties_solids", + "properties_vapor", + ] + + for port_str in port_lst: + assert hasattr(m.fs.unit, port_str) + port = getattr(m.fs.unit, port_str) + assert len(port.vars) == 3 + assert isinstance(port, Port) + for var_str in port_vars_lst: + assert hasattr(port, var_str) + var = getattr(port, var_str) + assert isinstance(var, Var) + + for n, eff in m.fs.unit.effects.items(): + for b in state_blks: + assert hasattr(eff.effect, b) + sb = getattr(eff.effect, b) + assert isinstance(sb, NaClStateBlock) + assert sb[0].temperature.ub == 1000 + + effect_params = [ + "approach_temperature_heat_exchanger", + "dimensionless_crystal_length", + ] + effect_vars = [ + "crystal_growth_rate", + "crystal_median_length", + "crystallization_yield", + "dens_mass_magma", + "dens_mass_slurry", + "diameter_crystallizer", + "height_crystallizer", + "height_slurry", + "magma_circulation_flow_vol", + "pressure_operating", + "product_volumetric_solids_fraction", + "relative_supersaturation", + "souders_brown_constant", + "t_res", + "temperature_operating", + "volume_suspension", + "work_mechanical", + ] + + effect_exprs = [ + "delta_temperature", + "eq_max_allowable_velocity", + "eq_minimum_height_diameter_ratio", + "eq_vapor_space_height", + ] + + effect_constr = [ + "eq_mass_balance_constraints", + "eq_solubility_massfrac_equality_constraint", + "eq_removal_balance", + "eq_vol_fraction_solids", + "eq_dens_magma", + "eq_operating_pressure_constraint", + "eq_relative_supersaturation", + "eq_enthalpy_balance", + "eq_p_con3", + "eq_T_con1", + "eq_T_con2", + "eq_T_con3", + "eq_minimum_hex_circulation_rate_constraint", + "eq_dens_mass_slurry", + "eq_residence_time", + "eq_suspension_volume", + "eq_vapor_head_diameter_constraint", + "eq_slurry_height_constraint", + "eq_crystallizer_height_constraint", + "eq_pure_water_mass_flow_rate", + "eq_vapor_energy_constraint", + "eq_p_con1", + "eq_p_con2", + "eq_p_con4", + "eq_p_con5", + ] + + assert hasattr(m.fs.unit, "recovery_vol_phase") + assert isinstance(m.fs.unit.recovery_vol_phase, Var) + + for n, eff in m.fs.unit.effects.items(): + for p in effect_params: + assert hasattr(eff.effect, p) + for v in effect_vars: + assert hasattr(eff.effect, v) + for e in effect_exprs: + assert hasattr(eff.effect, e) + for c in effect_constr: + assert hasattr(eff.effect, c) + assert hasattr(eff.effect, f"eq_delta_temperature_inlet_effect_{n}") + assert hasattr(eff.effect, f"eq_delta_temperature_outlet_effect_{n}") + assert hasattr(eff.effect, f"eq_heat_transfer_effect_{n}") + if n == 1: + assert number_variables(eff.effect) == 154 + assert number_total_constraints(eff.effect) == 128 + assert number_unused_variables(eff.effect) == 2 + assert hasattr(eff.effect, "heating_steam") + assert isinstance(eff.effect.heating_steam, WaterStateBlock) + assert eff.effect.heating_steam[0].temperature.ub == 1000 + assert hasattr(eff.effect, "eq_heating_steam_flow_rate") + if n != 1: + assert number_variables(eff.effect) == 148 + assert number_total_constraints(eff.effect) == 126 + assert number_unused_variables(eff.effect) == 0 + assert hasattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + + assert number_variables(m) == 609 + assert number_total_constraints(m) == 393 + assert number_unused_variables(m) == 43 + + assert_units_consistent(m) + + @pytest.mark.unit + def test_dof(self, MEC3_frame): + m = MEC3_frame + # With mass flow rates into each of the effects fixed, + # the model is over specified + assert degrees_of_freedom(m) == -2 + # Unfixing the mass flow rates for CV will result in 0 DOF + # (This is only done for testing purposes) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "H2O" + ].unfix() + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "NaCl" + ].unfix() + assert degrees_of_freedom(m) == 0 + for n, eff in m.fs.unit.effects.items(): + if n == 1: + assert degrees_of_freedom(eff.effect) == 0 + else: + assert degrees_of_freedom(eff.effect) == 3 + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "H2O" + ].fix() + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "NaCl" + ].fix() + assert degrees_of_freedom(m) == -2 + # Alternatively, one can .set_value() each of the mass flows for the individual effects, + # resulting in 4 DOF: + for n, eff in m.fs.unit.effects.items(): + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].unfix() + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].unfix() + assert degrees_of_freedom(m) == 4 + + @pytest.mark.component + def test_calculate_scaling(self, MEC3_frame): + m = MEC3_frame + + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "NaCl") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Vap", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Sol", "NaCl") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Vap", "H2O") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + + calculate_scaling_factors(m) + unscaled_var_list = list(unscaled_variables_generator(m)) + assert len(unscaled_var_list) == 0 + badly_scaled_var_list = list(badly_scaled_var_generator(m)) + assert len(badly_scaled_var_list) == 0 + + @pytest.mark.component + def test_initialize(self, MEC3_frame): + m = MEC3_frame + m.fs.unit.initialize() + c0 = value( + m.fs.unit.effects[1] + .effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + ) + htc = value(m.fs.unit.effects[1].effect.overall_heat_transfer_coefficient) + assert degrees_of_freedom(m) == 0 + for n, eff in m.fs.unit.effects.items(): + assert ( + not eff.effect.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .is_fixed() + ) + assert ( + not eff.effect.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .is_fixed() + ) + assert ( + pytest.approx( + value( + eff.effect.properties_in[0].conc_mass_phase_comp["Liq", "NaCl"] + ), + rel=1e-6, + ) + == c0 + ) + if n == 1: + assert degrees_of_freedom(eff.effect) == 2 + if n != 1: + assert ( + eff.effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + .is_fixed() + ) + linking_constr = getattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + assert linking_constr.active + assert eff.effect.overall_heat_transfer_coefficient.is_fixed() + assert value(eff.effect.overall_heat_transfer_coefficient) == htc + assert degrees_of_freedom(eff.effect) == 3 + + @pytest.mark.component + def test_solve(self, MEC3_frame): + + m = MEC3_frame + + results = solver.solve(m) + assert_optimal_termination(results) + + @pytest.mark.component + def test_conservation(self, MEC3_frame): + m = MEC3_frame + + comp_lst = ["NaCl", "H2O"] + phase_lst = ["Sol", "Liq", "Vap"] + + total_mass_flow_water_in = 0 + total_mass_flow_salt_in = 0 + total_mass_flow_water_out = 0 + + # Test mass balance for each effect + for _, e in m.fs.unit.effects.items(): + eff = e.effect + + phase_comp_list = [ + (p, j) + for j in comp_lst + for p in phase_lst + if (p, j) in eff.properties_in[0].phase_component_set + ] + flow_mass_in = sum( + eff.properties_in[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_out = sum( + eff.properties_out[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_solids = sum( + eff.properties_solids[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_vapor = sum( + eff.properties_vapor[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + + assert ( + abs( + value( + flow_mass_in + - flow_mass_out + - flow_mass_solids + - flow_mass_vapor + ) + ) + <= 1e-6 + ) + + total_mass_flow_water_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + ) + total_mass_flow_salt_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + total_mass_flow_water_out += value( + eff.properties_pure_water[0].flow_mass_phase_comp["Liq", "H2O"] + ) + + # Test control volume mass balance + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .value, + rel=1e-6, + ) + == total_mass_flow_salt_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_out[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_out + ) + + @pytest.mark.component + def test_solution(self, MEC3_frame): + m = MEC3_frame + + assert ( + pytest.approx(value(m.fs.unit.recovery_vol_phase["Liq"]), rel=1e-3) + == 0.5485 + ) + + unit_results_dict = { + 1: { + "product_volumetric_solids_fraction": 0.15774, + "temperature_operating": 359.4, + "pressure_operating": 45000.0, + "dens_mass_magma": 333.6, + "dens_mass_slurry": 1321.5, + "work_mechanical": {0.0: 13035.6}, + "diameter_crystallizer": 2.970, + "height_slurry": 2.377, + "height_crystallizer": 4.605, + "magma_circulation_flow_vol": 0.92570, + "relative_supersaturation": {"NaCl": 0.66123}, + "t_res": 1.02282118, + "volume_suspension": 16.4746, + "eq_vapor_space_height": 2.227, + "eq_minimum_height_diameter_ratio": 4.455, + "energy_flow_superheated_vapor": 11486.5, + "delta_temperature_in": {0.0: 61.67}, + "delta_temperature_out": {0.0: 128.0}, + "delta_temperature": {0.0: 90.8}, + "heat_exchanger_area": 1435.5, + }, + 2: { + "product_volumetric_solids_fraction": 0.156676, + "temperature_operating": 344.8, + "pressure_operating": 25000.0, + "dens_mass_magma": 331.3, + "dens_mass_slurry": 1324.8, + "work_mechanical": {0.0: 11486.5}, + "diameter_crystallizer": 3.234, + "height_slurry": 1.846, + "height_crystallizer": 4.85, + "magma_circulation_flow_vol": 0.81262, + "relative_supersaturation": {"NaCl": 0.66644}, + "t_res": 1.0228, + "volume_suspension": 15.1726, + "eq_minimum_height_diameter_ratio": 4.851, + "energy_flow_superheated_vapor": 10540.2, + "delta_temperature_in": {0.0: 14.61}, + "delta_temperature_out": {0.0: 58.77}, + "delta_temperature": {0.0: 31.59}, + "heat_exchanger_area": 3635.8, + }, + 3: { + "product_volumetric_solids_fraction": 0.15640, + "temperature_operating": 340.5, + "pressure_operating": 20800.0, + "dens_mass_magma": 330.8, + "dens_mass_slurry": 1325.96, + "work_mechanical": {0.0: 10540.2}, + "diameter_crystallizer": 3.2458, + "height_slurry": 1.704, + "height_crystallizer": 4.868, + "magma_circulation_flow_vol": 0.744761, + "relative_supersaturation": {"NaCl": 0.66785}, + "t_res": 1.022, + "volume_suspension": 14.10, + "energy_flow_superheated_vapor": 9791.82, + "delta_temperature_in": {0.0: 4.342}, + "delta_temperature_out": {0.0: 44.8}, + "delta_temperature": {0.0: 16.85}, + "heat_exchanger_area": 6254.11, + }, + } + + for n, d in unit_results_dict.items(): + eff = m.fs.unit.effects[n].effect + for v, r in d.items(): + effv = getattr(eff, v) + if effv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(effv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(effv), rel=1e-3) == r + + steam_results_dict = { + "flow_mass_phase_comp": {("Vap", "H2O"): 6.14896}, + "temperature": 421.1, + "pressure": 451325.0, + "dh_vap_mass": 2119982.2, + "pressure_sat": 451324.9, + } + + for v, r in steam_results_dict.items(): + sv = getattr(m.fs.unit.effects[1].effect.heating_steam[0], v) + if sv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(sv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(sv), rel=1e-3) == r + + @pytest.mark.component + def test_costing(self, MEC3_frame): + m = MEC3_frame + m.fs.costing = TreatmentCosting() + # m.fs.costing.base_currency = pyunits.USD_2018 + m.fs.unit.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method_arguments={"cost_type": "mass_basis"}, + ) + + m.fs.costing.nacl_recovered.cost.set_value(-0.024) + m.fs.costing.cost_process() + m.fs.costing.add_LCOW(m.fs.unit.total_flow_vol_in) + m.fs.costing.add_specific_energy_consumption( + m.fs.unit.total_flow_vol_in, name="SEC" + ) + results = solver.solve(m) + assert_optimal_termination(results) + + sys_costing_dict = { + "aggregate_capital_cost": 20645631.14, + "aggregate_flow_electricity": 46.05, + "aggregate_flow_NaCl_recovered": 4.1249, + "aggregate_flow_steam": 2.6482, + "aggregate_flow_costs": { + "electricity": 33176.93, + "NaCl_recovered": -3710055.74, + "steam": 392435.52, + }, + "aggregate_direct_capital_cost": 10322815.57, + "total_capital_cost": 20645631.14, + "total_operating_cost": -2665074.35, + "maintenance_labor_chemical_operating_cost": 619368.93, + "total_fixed_operating_cost": 619368.93, + "total_variable_operating_cost": -3284443.28, + "total_annualized_cost": -353673.11, + "LCOW": -0.443848, + "SEC": 0.506672, + } + + for v, r in sys_costing_dict.items(): + cv = getattr(m.fs.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + eff_costing_dict = { + "capital_cost": 20645631.14, + "direct_capital_cost": 10322815.57, + "capital_cost_crystallizer_effect_1": 3216728.22, + "capital_cost_heat_exchanger_effect_1": 1462729.88, + "capital_cost_effect_1": 2339729.05, + "capital_cost_crystallizer_effect_2": 3068214.17, + "capital_cost_heat_exchanger_effect_2": 3667644.07, + "capital_cost_effect_2": 3367929.12, + "capital_cost_crystallizer_effect_3": 2949086.04, + "capital_cost_heat_exchanger_effect_3": 6281228.74, + "capital_cost_effect_3": 4615157.39, + } + + for v, r in eff_costing_dict.items(): + cv = getattr(m.fs.unit.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + +class TestMultiEffectCrystallizer_4Effects: + @pytest.fixture(scope="class") + def MEC4_frame(self): + """ + Test crystallizer with 4 effects + """ + m = build_mec4() + return m + + @pytest.mark.unit + def test_config(self, MEC4_frame): + m = MEC4_frame + assert len(m.fs.unit.config) == 6 + + assert not m.fs.unit.config.dynamic + assert not m.fs.unit.config.has_holdup + assert m.fs.unit.config.property_package is m.fs.properties + assert m.fs.unit.config.property_package_vapor is m.fs.vapor_properties + assert m.fs.unit.config.number_effects == 4 + assert m.fs.unit.config.number_effects == len(m.fs.unit.effects) + + for _, eff in m.fs.unit.effects.items(): + assert isinstance(eff.effect, CrystallizerEffect) + assert not eff.effect.config.standalone + + @pytest.mark.unit + def test_build(self, MEC4_frame): + m = MEC4_frame + + # test ports and variables + port_lst = ["inlet", "outlet", "solids", "vapor", "steam"] + port_vars_lst = ["flow_mass_phase_comp", "pressure", "temperature"] + state_blks = [ + "properties_in", + "properties_out", + "properties_pure_water", + "properties_solids", + "properties_vapor", + ] + + for port_str in port_lst: + assert hasattr(m.fs.unit, port_str) + port = getattr(m.fs.unit, port_str) + assert len(port.vars) == 3 + assert isinstance(port, Port) + for var_str in port_vars_lst: + assert hasattr(port, var_str) + var = getattr(port, var_str) + assert isinstance(var, Var) + + for n, eff in m.fs.unit.effects.items(): + for b in state_blks: + assert hasattr(eff.effect, b) + sb = getattr(eff.effect, b) + assert isinstance(sb, NaClStateBlock) + assert sb[0].temperature.ub == 1000 + + effect_params = [ + "approach_temperature_heat_exchanger", + "dimensionless_crystal_length", + ] + effect_vars = [ + "crystal_growth_rate", + "crystal_median_length", + "crystallization_yield", + "dens_mass_magma", + "dens_mass_slurry", + "diameter_crystallizer", + "height_crystallizer", + "height_slurry", + "magma_circulation_flow_vol", + "pressure_operating", + "product_volumetric_solids_fraction", + "relative_supersaturation", + "souders_brown_constant", + "t_res", + "temperature_operating", + "volume_suspension", + "work_mechanical", + ] + + effect_exprs = [ + "delta_temperature", + "eq_max_allowable_velocity", + "eq_minimum_height_diameter_ratio", + "eq_vapor_space_height", + ] + + effect_constr = [ + "eq_mass_balance_constraints", + "eq_solubility_massfrac_equality_constraint", + "eq_removal_balance", + "eq_vol_fraction_solids", + "eq_dens_magma", + "eq_operating_pressure_constraint", + "eq_relative_supersaturation", + "eq_enthalpy_balance", + "eq_p_con3", + "eq_T_con1", + "eq_T_con2", + "eq_T_con3", + "eq_minimum_hex_circulation_rate_constraint", + "eq_dens_mass_slurry", + "eq_residence_time", + "eq_suspension_volume", + "eq_vapor_head_diameter_constraint", + "eq_slurry_height_constraint", + "eq_crystallizer_height_constraint", + "eq_pure_water_mass_flow_rate", + "eq_vapor_energy_constraint", + "eq_p_con1", + "eq_p_con2", + "eq_p_con4", + "eq_p_con5", + ] + + assert hasattr(m.fs.unit, "recovery_vol_phase") + assert isinstance(m.fs.unit.recovery_vol_phase, Var) + + for n, eff in m.fs.unit.effects.items(): + for p in effect_params: + assert hasattr(eff.effect, p) + for v in effect_vars: + assert hasattr(eff.effect, v) + for e in effect_exprs: + assert hasattr(eff.effect, e) + for c in effect_constr: + assert hasattr(eff.effect, c) + assert hasattr(eff.effect, f"eq_delta_temperature_inlet_effect_{n}") + assert hasattr(eff.effect, f"eq_delta_temperature_outlet_effect_{n}") + assert hasattr(eff.effect, f"eq_heat_transfer_effect_{n}") + if n == 1: + assert number_variables(eff.effect) == 154 + assert number_total_constraints(eff.effect) == 128 + assert number_unused_variables(eff.effect) == 2 + assert hasattr(eff.effect, "heating_steam") + assert isinstance(eff.effect.heating_steam, WaterStateBlock) + assert eff.effect.heating_steam[0].temperature.ub == 1000 + assert hasattr(eff.effect, "eq_heating_steam_flow_rate") + if n != 1: + assert number_variables(eff.effect) == 148 + assert number_total_constraints(eff.effect) == 126 + assert number_unused_variables(eff.effect) == 0 + assert hasattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + + assert number_variables(m) == 757 + assert number_total_constraints(m) == 519 + assert number_unused_variables(m) == 43 + + assert_units_consistent(m) + + @pytest.mark.unit + def test_dof(self, MEC4_frame): + + m = MEC4_frame + # With mass flow rates into each of the effects fixed, + # and the control_volume flows fixed, the model is over specified + assert degrees_of_freedom(m) == -2 + # Unfixing the mass flow rates for CV will result in 0 DOF + # (This is only done for testing purposes) + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "H2O" + ].unfix() + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "NaCl" + ].unfix() + assert degrees_of_freedom(m) == 0 + for n, eff in m.fs.unit.effects.items(): + if n == 1: + assert degrees_of_freedom(eff.effect) == 0 + else: + assert degrees_of_freedom(eff.effect) == 3 + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "H2O" + ].fix() + m.fs.unit.control_volume.properties_in[0].flow_mass_phase_comp[ + "Liq", "NaCl" + ].fix() + assert degrees_of_freedom(m) == -2 + # Alternatively, one can .set_value() each of the mass flows for the individual effects, + # resulting in 6 DOF: + for n, eff in m.fs.unit.effects.items(): + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].unfix() + eff.effect.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"].unfix() + assert degrees_of_freedom(m) == 6 + + @pytest.mark.unit + def test_calculate_scaling(self, MEC4_frame): + m = MEC4_frame + + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Liq", "NaCl") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Vap", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Sol", "NaCl") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1e-1, index=("Vap", "H2O") + ) + m.fs.vapor_properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + + calculate_scaling_factors(m) + unscaled_var_list = list(unscaled_variables_generator(m)) + assert len(unscaled_var_list) == 0 + badly_scaled_var_list = list(badly_scaled_var_generator(m)) + assert len(badly_scaled_var_list) == 0 + + @pytest.mark.component + def test_initialize(self, MEC4_frame): + m = MEC4_frame + m.fs.unit.initialize() + c0 = value( + m.fs.unit.effects[1] + .effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + ) + htc = value(m.fs.unit.effects[1].effect.overall_heat_transfer_coefficient) + assert degrees_of_freedom(m) == 0 + for n, eff in m.fs.unit.effects.items(): + assert ( + not eff.effect.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .is_fixed() + ) + assert ( + not eff.effect.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .is_fixed() + ) + assert ( + pytest.approx( + value( + eff.effect.properties_in[0].conc_mass_phase_comp["Liq", "NaCl"] + ), + rel=1e-6, + ) + == c0 + ) + if n == 1: + assert degrees_of_freedom(eff.effect) == 2 + if n != 1: + assert ( + eff.effect.properties_in[0] + .conc_mass_phase_comp["Liq", "NaCl"] + .is_fixed() + ) + linking_constr = getattr( + eff.effect, f"eq_energy_for_effect_{n}_from_effect_{n - 1}" + ) + assert linking_constr.active + assert eff.effect.overall_heat_transfer_coefficient.is_fixed() + assert value(eff.effect.overall_heat_transfer_coefficient) == htc + assert degrees_of_freedom(eff.effect) == 3 + + @pytest.mark.component + def test_solve(self, MEC4_frame): + + m = MEC4_frame + + results = solver.solve(m) + assert_optimal_termination(results) + + @pytest.mark.component + def test_conservation(self, MEC4_frame): + + m = MEC4_frame + + comp_lst = ["NaCl", "H2O"] + phase_lst = ["Sol", "Liq", "Vap"] + + total_mass_flow_water_in = 0 + total_mass_flow_salt_in = 0 + total_mass_flow_water_out = 0 + + for _, e in m.fs.unit.effects.items(): + eff = e.effect + + phase_comp_list = [ + (p, j) + for j in comp_lst + for p in phase_lst + if (p, j) in eff.properties_in[0].phase_component_set + ] + flow_mass_in = sum( + eff.properties_in[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_out = sum( + eff.properties_out[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_solids = sum( + eff.properties_solids[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_vapor = sum( + eff.properties_vapor[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + + assert ( + abs( + value( + flow_mass_in + - flow_mass_out + - flow_mass_solids + - flow_mass_vapor + ) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + flow_mass_in * eff.properties_in[0].enth_mass_phase["Liq"] + - flow_mass_out * eff.properties_out[0].enth_mass_phase["Liq"] + - flow_mass_vapor + * eff.properties_vapor[0].enth_mass_solvent["Vap"] + - flow_mass_solids + * eff.properties_solids[0].enth_mass_solute["Sol"] + - flow_mass_solids + * eff.properties_solids[0].dh_crystallization_mass_comp["NaCl"] + + eff.work_mechanical[0] + ) + ) + <= 1e-2 + ) + + total_mass_flow_water_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + ) + total_mass_flow_salt_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + total_mass_flow_water_out += value( + eff.properties_pure_water[0].flow_mass_phase_comp["Liq", "H2O"] + ) + + # Test control volume mass balance + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .value, + rel=1e-6, + ) + == total_mass_flow_salt_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_out[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_out + ) + + @pytest.mark.component + def test_solution(self, MEC4_frame): + m = MEC4_frame + + assert ( + pytest.approx(value(m.fs.unit.recovery_vol_phase["Liq"]), rel=1e-3) + == 0.734452 + ) + + unit_results_dict = { + 1: { + "crystallization_yield": {"NaCl": 0.5}, + "product_volumetric_solids_fraction": 0.132962, + "temperature_operating": 359.48, + "pressure_operating": 45000.0, + "dens_mass_magma": 281.21, + "dens_mass_slurry": 1298.22, + "work_mechanical": {0.0: 1915.44}, + "diameter_crystallizer": 1.1448, + "height_slurry": 1.073, + "height_crystallizer": 1.9316, + "magma_circulation_flow_vol": 0.134172, + "relative_supersaturation": {"NaCl": 0.567033}, + "t_res": 1.0228, + "volume_suspension": 1.1045, + "eq_max_allowable_velocity": 2.6304, + "eq_vapor_space_height": 0.858635, + "eq_minimum_height_diameter_ratio": 1.7172, + "energy_flow_superheated_vapor": 1706.7, + "delta_temperature_in": {0.0: 57.39}, + "delta_temperature_out": {0.0: 123.72}, + "delta_temperature": {0.0: 86.31}, + "heat_exchanger_area": 221.91, + }, + 2: { + "product_volumetric_solids_fraction": 0.132109, + "temperature_operating": 344.86, + "pressure_operating": 25000.0, + "dens_mass_magma": 279.41, + "dens_mass_slurry": 1301.84, + "work_mechanical": {0.0: 1706.7}, + "diameter_crystallizer": 1.2481, + "height_slurry": 0.828632, + "height_crystallizer": 1.8721, + "magma_circulation_flow_vol": 0.119101, + "relative_supersaturation": {"NaCl": 0.571246}, + "t_res": 1.0228, + "volume_suspension": 1.0138, + "eq_max_allowable_velocity": 3.4641, + "eq_vapor_space_height": 0.936079, + "eq_minimum_height_diameter_ratio": 1.8721, + "energy_flow_superheated_vapor": 1569.57, + "delta_temperature_in": {0.0: 14.61}, + "delta_temperature_out": {0.0: 58.77}, + "delta_temperature": {0.0: 31.59}, + "heat_exchanger_area": 540.22, + }, + 3: { + "product_volumetric_solids_fraction": 0.131922, + "temperature_operating": 340.52, + "pressure_operating": 20800.0, + "dens_mass_magma": 279.01, + "dens_mass_slurry": 1303.05, + "work_mechanical": {0.0: 1569.57}, + "diameter_crystallizer": 1.2521, + "height_slurry": 0.763702, + "height_crystallizer": 1.8782, + "magma_circulation_flow_vol": 0.109397, + "relative_supersaturation": {"NaCl": 0.572391}, + "t_res": 1.0228, + "volume_suspension": 0.940428, + "eq_max_allowable_velocity": 3.7763, + "eq_vapor_space_height": 0.939111, + "eq_minimum_height_diameter_ratio": 1.8782, + "energy_flow_superheated_vapor": 1457.19, + "delta_temperature_in": {0.0: 4.3421}, + "delta_temperature_out": {0.0: 44.83}, + "delta_temperature": {0.0: 16.85}, + "heat_exchanger_area": 931.31, + }, + 4: { + "product_volumetric_solids_fraction": 0.131442, + "temperature_operating": 323.22, + "pressure_operating": 9500.0, + "dens_mass_magma": 278.0, + "dens_mass_slurry": 1308.49, + "work_mechanical": {0.0: 1457.19}, + "diameter_crystallizer": 1.4623, + "height_slurry": 0.537417, + "height_crystallizer": 2.1935, + "magma_circulation_flow_vol": 0.101015, + "relative_supersaturation": {"NaCl": 0.576466}, + "t_res": 1.0228, + "volume_suspension": 0.902678, + "eq_max_allowable_velocity": 5.4596, + "eq_vapor_space_height": 1.0967, + "eq_minimum_height_diameter_ratio": 2.1935, + "energy_flow_superheated_vapor": 1405.09, + "delta_temperature_in": {0.0: 17.29}, + "delta_temperature_out": {0.0: 40.69}, + "delta_temperature": {0.0: 27.32}, + "heat_exchanger_area": 533.18, + }, + } + + for n, d in unit_results_dict.items(): + eff = m.fs.unit.effects[n].effect + for v, r in d.items(): + effv = getattr(eff, v) + if effv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(effv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(effv), rel=1e-3) == r + + steam_results_dict = { + "flow_mass_phase_comp": {("Liq", "H2O"): 0.0, ("Vap", "H2O"): 0.89795}, + "temperature": 416.8, + "pressure": 401325.0, + "dh_vap_mass": 2133119.0, + "pressure_sat": 401325.0, + } + + for v, r in steam_results_dict.items(): + sv = getattr(m.fs.unit.effects[1].effect.heating_steam[0], v) + if sv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(sv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(sv), rel=1e-3) == r + + @pytest.mark.component + def test_costing(self, MEC4_frame): + m = MEC4_frame + m.fs.costing = TreatmentCosting() + m.fs.unit.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + ) + + m.fs.costing.nacl_recovered.cost.set_value(-0.024) + m.fs.costing.cost_process() + m.fs.costing.add_LCOW(m.fs.unit.total_flow_vol_in) + m.fs.costing.add_specific_energy_consumption( + m.fs.unit.total_flow_vol_in, name="SEC" + ) + results = solver.solve(m) + assert_optimal_termination(results) + + sys_costing_dict = { + "aggregate_capital_cost": 4934127.42, + "aggregate_flow_electricity": 8.4612, + "aggregate_flow_NaCl_recovered": 0.299999, + "aggregate_flow_steam": 0.4304, + "aggregate_flow_costs": { + "electricity": 6095.09, + "NaCl_recovered": -269822.09, + "steam": 63793.11, + }, + "aggregate_direct_capital_cost": 2467063.71, + "total_capital_cost": 4934127.42, + "total_operating_cost": -51910.07, + "maintenance_labor_chemical_operating_cost": 148023.82, + "total_fixed_operating_cost": 148023.82, + "total_variable_operating_cost": -199933.89, + "total_annualized_cost": 500494.84, + "LCOW": 4.3957, + "SEC": 0.65144, + } + + for v, r in sys_costing_dict.items(): + cv = getattr(m.fs.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + eff_costing_dict = { + "capital_cost": 4934127.42, + "cost_factor": 2.0, + "direct_capital_cost": 2467063.71, + "capital_cost_crystallizer_effect_1": 701221.05, + "capital_cost_heat_exchanger_effect_1": 234213.59, + "capital_cost_effect_1": 467717.32, + "capital_cost_crystallizer_effect_2": 667484.84, + "capital_cost_heat_exchanger_effect_2": 558948.69, + "capital_cost_effect_2": 613216.76, + "capital_cost_crystallizer_effect_3": 640779.56, + "capital_cost_heat_exchanger_effect_3": 954745.08, + "capital_cost_effect_3": 797762.32, + "capital_cost_crystallizer_effect_4": 624930.83, + "capital_cost_heat_exchanger_effect_4": 551803.75, + "capital_cost_effect_4": 588367.29, + } + + for v, r in eff_costing_dict.items(): + cv = getattr(m.fs.unit.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + @pytest.mark.component + def test_costing_by_volume(self): + + m = build_mec4() + + for _, eff in m.fs.unit.effects.items(): + eff.effect.crystal_median_length.fix(0.6e-3) + eff.effect.crystal_growth_rate.fix(5e-9) + + self.test_calculate_scaling(m) + + m.fs.unit.initialize() + + assert degrees_of_freedom(m) == 0 + + results = solver.solve(m) + assert_optimal_termination(results) + + m.fs.costing = TreatmentCosting() + m.fs.unit.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method_arguments={"cost_type": "volume_basis"}, + ) + + m.fs.costing.nacl_recovered.cost.set_value(-0.024) + m.fs.costing.cost_process() + m.fs.costing.add_LCOW(m.fs.unit.total_flow_vol_in) + m.fs.costing.add_specific_energy_consumption( + m.fs.unit.total_flow_vol_in, name="SEC" + ) + results = solver.solve(m) + + assert_optimal_termination(results) + + sys_costing_dict = { + "aggregate_capital_cost": 5078402.93, + "aggregate_flow_electricity": 8.4612, + "aggregate_flow_NaCl_recovered": 0.299999, + "aggregate_flow_steam": 0.430492, + "aggregate_flow_costs": { + "electricity": 6095.09, + "NaCl_recovered": -269822.09, + "steam": 63793.11, + }, + "aggregate_direct_capital_cost": 2539201.46, + "total_capital_cost": 5078402.93, + "total_operating_cost": -47581.81, + "maintenance_labor_chemical_operating_cost": 152352.08, + "total_fixed_operating_cost": 152352.08, + "total_variable_operating_cost": -199933.89, + "total_annualized_cost": 520975.61, + "LCOW": 4.5756, + "SEC": 0.65144, + } + + for v, r in sys_costing_dict.items(): + cv = getattr(m.fs.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + eff_costing_dict = { + "capital_cost": 5078402.93, + "cost_factor": 2.0, + "direct_capital_cost": 2539201.46, + "capital_cost_crystallizer_effect_1": 715324.02, + "capital_cost_heat_exchanger_effect_1": 234213.59, + "capital_cost_effect_1": 474768.8, + "capital_cost_crystallizer_effect_2": 697956.34, + "capital_cost_heat_exchanger_effect_2": 558948.69, + "capital_cost_effect_2": 628452.52, + "capital_cost_crystallizer_effect_3": 676895.98, + "capital_cost_heat_exchanger_effect_3": 954745.08, + "capital_cost_effect_3": 815820.53, + "capital_cost_crystallizer_effect_4": 688515.45, + "capital_cost_heat_exchanger_effect_4": 551803.75, + "capital_cost_effect_4": 620159.6, + } + + for v, r in eff_costing_dict.items(): + cv = getattr(m.fs.unit.costing, v) + if cv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(cv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(cv), rel=1e-3) == r + + @pytest.mark.component + def test_optimization(self, MEC4_frame): + m = MEC4_frame + + # Optimization scenario + # Looking for the best operating pressures within a typical range + + for n, eff in m.fs.unit.effects.items(): + eff.effect.pressure_operating.unfix() + if n == m.fs.unit.first_effect: + eff.effect.pressure_operating.setub(1.2e5) + if n == m.fs.unit.last_effect: + eff.effect.pressure_operating.setlb(0.02e5) + + eff_idx = [n for n in m.fs.unit.Effects if n != m.fs.unit.first_effect] + + @m.fs.unit.Constraint(eff_idx, doc="Pressure decreasing across effects") + def eq_pressure_decreasing_across_eff(b, n): + return ( + b.effects[n].effect.pressure_operating + <= b.effects[n - 1].effect.pressure_operating + ) + + m.fs.unit.temperature_diff_typical = Var( + initialize=12, + bounds=(0, None), + units=pyunits.degK, + doc="Typical temperature difference limit between effects in industry", + ) + + m.fs.unit.temperature_diff_typical.fix() + + @m.fs.unit.Constraint( + eff_idx, + doc="Temperature difference limit (based on industrial convention)", + ) + def eq_temperature_difference_limit(b, n): + return ( + b.effects[n].effect.temperature_operating + >= b.effects[n - 1].effect.temperature_operating + - b.temperature_diff_typical + ) + + m.fs.objective = Objective(expr=m.fs.costing.LCOW) + + opt_results = solver.solve(m) + assert_optimal_termination(opt_results) + + @pytest.mark.component + def test_optimization_solution(self, MEC4_frame): + m = MEC4_frame + + # Check the optimized LCOW + assert pytest.approx(value(m.fs.costing.LCOW), rel=1e-3) == 3.865 + + assert ( + pytest.approx(value(m.fs.unit.recovery_vol_phase["Liq"]), rel=1e-3) + == 0.7544 + ) + + unit_results_dict = { + 1: { + "product_volumetric_solids_fraction": 0.135236, + "temperature_operating": 387.01, + "pressure_operating": 120000.0, + "work_mechanical": {0.0: 2163.65}, + "heat_exchanger_area": 329.25, + }, + 2: { + "product_volumetric_solids_fraction": 0.134156, + "temperature_operating": 375.01, + "pressure_operating": 79812.1, + "work_mechanical": {0.0: 1818.52}, + "heat_exchanger_area": 495.94, + }, + 3: { + "product_volumetric_solids_fraction": 0.13323, + "temperature_operating": 363.01, + "pressure_operating": 51490.6, + "work_mechanical": {0.0: 1567.57}, + "heat_exchanger_area": 467.40, + }, + 4: { + "product_volumetric_solids_fraction": 0.132473, + "temperature_operating": 351.01, + "pressure_operating": 32193.7, + "work_mechanical": {0.0: 1386.22}, + "heat_exchanger_area": 458.50, + }, + } + + for n, d in unit_results_dict.items(): + eff = m.fs.unit.effects[n].effect + for v, r in d.items(): + effv = getattr(eff, v) + if effv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(effv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(effv), rel=1e-3) == r + + steam_results_dict = { + "flow_mass_phase_comp": {("Liq", "H2O"): 0.0, ("Vap", "H2O"): 1.014}, + "temperature": 416.8, + "pressure": 401325.0, + "dh_vap_mass": 2133119.0, + "pressure_sat": 401325.0, + } + + for v, r in steam_results_dict.items(): + sv = getattr(m.fs.unit.effects[1].effect.heating_steam[0], v) + if sv.is_indexed(): + for i, s in r.items(): + assert pytest.approx(value(sv[i]), rel=1e-3) == s + else: + assert pytest.approx(value(sv), rel=1e-3) == r + + @pytest.mark.component + def test_optimization_conservation(self, MEC4_frame): + + m = MEC4_frame + + comp_lst = ["NaCl", "H2O"] + phase_lst = ["Sol", "Liq", "Vap"] + + total_mass_flow_water_in = 0 + total_mass_flow_salt_in = 0 + total_mass_flow_water_out = 0 + + for _, e in m.fs.unit.effects.items(): + eff = e.effect + + phase_comp_list = [ + (p, j) + for j in comp_lst + for p in phase_lst + if (p, j) in eff.properties_in[0].phase_component_set + ] + flow_mass_in = sum( + eff.properties_in[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_out = sum( + eff.properties_out[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_solids = sum( + eff.properties_solids[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + flow_mass_vapor = sum( + eff.properties_vapor[0].flow_mass_phase_comp[p, j] + for p in phase_lst + for j in comp_lst + if (p, j) in phase_comp_list + ) + + assert ( + abs( + value( + flow_mass_in + - flow_mass_out + - flow_mass_solids + - flow_mass_vapor + ) + ) + <= 1e-6 + ) + + assert ( + abs( + value( + flow_mass_in * eff.properties_in[0].enth_mass_phase["Liq"] + - flow_mass_out * eff.properties_out[0].enth_mass_phase["Liq"] + - flow_mass_vapor + * eff.properties_vapor[0].enth_mass_solvent["Vap"] + - flow_mass_solids + * eff.properties_solids[0].enth_mass_solute["Sol"] + - flow_mass_solids + * eff.properties_solids[0].dh_crystallization_mass_comp["NaCl"] + + eff.work_mechanical[0] + ) + ) + <= 1e-2 + ) + + total_mass_flow_water_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + ) + total_mass_flow_salt_in += value( + eff.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + total_mass_flow_water_out += value( + eff.properties_pure_water[0].flow_mass_phase_comp["Liq", "H2O"] + ) + + # Test control volume mass balance + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_in[0] + .flow_mass_phase_comp["Liq", "NaCl"] + .value, + rel=1e-6, + ) + == total_mass_flow_salt_in + ) + assert ( + pytest.approx( + m.fs.unit.control_volume.properties_out[0] + .flow_mass_phase_comp["Liq", "H2O"] + .value, + rel=1e-6, + ) + == total_mass_flow_water_out + )