diff --git a/watertap/costing/unit_models/crystallizer.py b/watertap/costing/unit_models/crystallizer.py index 4c5e6939b9..66e9ef448d 100644 --- a/watertap/costing/unit_models/crystallizer.py +++ b/watertap/costing/unit_models/crystallizer.py @@ -93,7 +93,7 @@ def build_crystallizer_cost_param_block(blk): ) costing = blk.parent_block() - costing.register_flow_type("steam", blk.steam_cost) + #costing.register_flow_type("steam", blk.steam_cost) costing.register_flow_type("NaCl", blk.NaCl_recovery_value) @@ -120,27 +120,6 @@ def cost_crystallizer(blk, cost_type=CrystallizerCostType.default): def _cost_crystallizer_flows(blk): - blk.costing_package.cost_flow( - pyo.units.convert( - ( - blk.unit_model.magma_circulation_flow_vol - * blk.unit_model.dens_mass_slurry - * Constants.acceleration_gravity - * blk.costing_package.crystallizer.pump_head_height - / blk.costing_package.crystallizer.efficiency_pump - ), - to_units=pyo.units.kW, - ), - "electricity", - ) - - blk.costing_package.cost_flow( - pyo.units.convert( - (blk.unit_model.work_mechanical[0] / _compute_steam_properties(blk)), - to_units=pyo.units.m**3 / pyo.units.s, - ), - "steam", - ) blk.costing_package.cost_flow( blk.unit_model.solids.flow_mass_phase_comp[0, "Sol", "NaCl"], @@ -215,72 +194,3 @@ def cost_crystallizer_by_volume(blk): ) ) _cost_crystallizer_flows(blk) - - -def _compute_steam_properties(blk): - """ - Function for computing saturated steam properties for thermal heating estimation. - - Args: - pressure_sat: Steam gauge pressure in bar - - Out: - Steam thermal capacity (latent heat of condensation * density) in kJ/m3 - """ - pressure_sat = blk.costing_package.crystallizer.steam_pressure - # 1. Compute saturation temperature of steam: computed from El-Dessouky expression - tsat_constants = [ - 42.6776 * pyo.units.K, - -3892.7 * pyo.units.K, - 1000 * pyo.units.kPa, - -9.48654 * pyo.units.dimensionless, - ] - psat = ( - pyo.units.convert(pressure_sat, to_units=pyo.units.kPa) - + 101.325 * pyo.units.kPa - ) - temperature_sat = tsat_constants[0] + tsat_constants[1] / ( - pyo.log(psat / tsat_constants[2]) + tsat_constants[3] - ) - - # 2. Compute latent heat of condensation/vaporization: computed from Sharqawy expression - t = temperature_sat - 273.15 * pyo.units.K - enth_mass_units = pyo.units.J / pyo.units.kg - t_inv_units = pyo.units.K**-1 - dh_constants = [ - 2.501e6 * enth_mass_units, - -2.369e3 * enth_mass_units * t_inv_units**1, - 2.678e-1 * enth_mass_units * t_inv_units**2, - -8.103e-3 * enth_mass_units * t_inv_units**3, - -2.079e-5 * enth_mass_units * t_inv_units**4, - ] - dh_vap = ( - dh_constants[0] - + dh_constants[1] * t - + dh_constants[2] * t**2 - + dh_constants[3] * t**3 - + dh_constants[4] * t**4 - ) - dh_vap = pyo.units.convert(dh_vap, to_units=pyo.units.kJ / pyo.units.kg) - - # 3. Compute specific volume: computed from Affandi expression (Eq 5) - t_critical = 647.096 * pyo.units.K - t_red = temperature_sat / t_critical # Reduced temperature - sp_vol_constants = [ - -7.75883 * pyo.units.dimensionless, - 3.23753 * pyo.units.dimensionless, - 2.05755 * pyo.units.dimensionless, - -0.06052 * pyo.units.dimensionless, - 0.00529 * pyo.units.dimensionless, - ] - log_sp_vol = ( - sp_vol_constants[0] - + sp_vol_constants[1] * (pyo.log(1 / t_red)) ** 0.4 - + sp_vol_constants[2] / (t_red**2) - + sp_vol_constants[3] / (t_red**4) - + sp_vol_constants[4] / (t_red**5) - ) - sp_vol = pyo.exp(log_sp_vol) * pyo.units.m**3 / pyo.units.kg - - # 4. Return specific energy: density * latent heat - return dh_vap / sp_vol diff --git a/watertap/costing/unit_models/heat_exchanger.py b/watertap/costing/unit_models/heat_exchanger.py index edecf79636..da13180920 100644 --- a/watertap/costing/unit_models/heat_exchanger.py +++ b/watertap/costing/unit_models/heat_exchanger.py @@ -34,7 +34,8 @@ def build_heat_exchanger_cost_param_block(blk): doc="steam cost per kg", ) - blk.parent_block().register_flow_type("steam", blk.steam_cost) + costing = blk.parent_block() + costing.register_flow_type("steam", blk.steam_cost) @register_costing_parameter_block( diff --git a/watertap/flowsheets/crystallization/Crystallizer_MVR.py b/watertap/flowsheets/crystallization/Crystallizer_MVR.py new file mode 100644 index 0000000000..f40a26a072 --- /dev/null +++ b/watertap/flowsheets/crystallization/Crystallizer_MVR.py @@ -0,0 +1,444 @@ +################################################################################# +# 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 pyomo.environ import ( + ConcreteModel, + TerminationCondition, +) +from pyomo.environ import ( + ConcreteModel, + value, + Constraint, + Objective, + Var, + NonNegativeReals, + TransformationFactory, + units as pyunits, + check_optimal_termination, +) +from pyomo.network import Arc +from idaes.core import FlowsheetBlock +from watertap.core.solvers import get_solver +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.initialization import ( + propagate_state, +) +from watertap.unit_models.pressure_changer import Pump +from watertap.costing.unit_models.pump import ( + cost_pump, +) +from pyomo.util.check_units import assert_units_consistent + +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom + +import idaes.core.util.scaling as iscale +import idaes.logger as idaeslog +from watertap.core.solvers import get_solver +from idaes.core import UnitModelCostingBlock + +from watertap.property_models.unit_specific import cryst_prop_pack as props +from watertap.unit_models.crystallizer import Crystallization +from watertap.costing import WaterTAPCosting, CrystallizerCostType + +from io import StringIO +from pyomo.util.infeasible import ( + log_infeasible_constraints, +) +from watertap.unit_models.steam_heater_0D import SteamHeater0D, Mode +from idaes.models.unit_models import Heater, Separator, Mixer, Product, Feed +from idaes.models.unit_models.mixer import MomentumMixingType +from watertap.core.solvers import get_solver +from watertap.unit_models.mvc.components.lmtd_chen_callback import ( + delta_temperature_chen_callback, +) +from idaes.models.unit_models.heat_exchanger import ( + HeatExchangerFlowPattern, +) +from watertap.unit_models.mvc.components import Compressor +from watertap.property_models.unit_specific import cryst_prop_pack as props +import watertap.property_models.water_prop_pack as props_w +import watertap.property_models.NaCl_T_dep_prop_pack as props_nacl +from pyomo.common.log import LoggingIntercept +from idaes.models.unit_models.translator import Translator +from idaes.models.unit_models.separator import SplittingType +import logging +#from watertap.core.util.model_debug_mode import activate +#activate() + + + +__author__ = "Elmira Shamlou" + + +def main(): + solver = get_solver() + m = build() + set_operating_conditions(m) + initialize_system(m, solver=solver) + + optimize_set_up(m) + + solve(m, solver=solver) + + display_system(m) + + + return m + + +def build(): + + # flowsheet set up + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + # properties + m.fs.properties_vapor = props_w.WaterParameterBlock() + m.fs.properties = props.NaClParameterBlock() + m.fs.properties_nacl = props_nacl.NaClParameterBlock() + + # Control volume flow blocks + m.fs.feed = Feed(property_package= m.fs.properties_nacl) + m.fs.distillate = Product(property_package= m.fs.properties_vapor) + + + # unit models: steam heater, mixer, pump, crystalizer, compressor, separator + + m.fs.crystallizer = Crystallization(property_package=m.fs.properties) + m.fs.heater = SteamHeater0D( + hot_side_name="hot", + cold_side_name="cold", + hot={ + "property_package": m.fs.properties_vapor, + }, + cold={ + "property_package": m.fs.properties_nacl, + }, + delta_temperature_callback=delta_temperature_chen_callback, + flow_pattern=HeatExchangerFlowPattern.countercurrent, + mode=Mode.CONDENSER, + ) + + m.fs.mixer = Mixer( + property_package=m.fs.properties_nacl, + momentum_mixing_type=MomentumMixingType.equality, + inlet_list=["feed", "recycle"], + ) + m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + m.fs.compressor = Compressor(property_package=m.fs.properties_vapor) + + m.fs.pump = Pump(property_package=m.fs.properties_nacl) + + m.fs.separator = Separator( + property_package=m.fs.properties_nacl, + outlet_list=["recycle", "purge"], + split_basis=SplittingType.totalFlow, + ) + + m.fs.tb_vapor = Translator( + inlet_property_package=m.fs.properties, + outlet_property_package=m.fs.properties_vapor, + ) + + @m.fs.tb_vapor.Constraint() + def eq_flow_mass_comp(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] + ) + @m.fs.tb_vapor.Constraint() + def eq_flow_mass_comp_vapor(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Vap", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Vap", "H2O"] + ) + + @m.fs.tb_vapor.Constraint() + def eq_temperature(blk): + return blk.properties_in[0].temperature == blk.properties_out[0].temperature + + @m.fs.tb_vapor.Constraint() + def eq_pressure(blk): + return blk.properties_in[0].pressure == blk.properties_out[0].pressure + + m.fs.tb_heater_to_cryst = Translator( + inlet_property_package=m.fs.properties_nacl, + outlet_property_package=m.fs.properties, + ) + + @m.fs.tb_heater_to_cryst.Constraint() + def eq_flow_mass_comp(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] + ) + @m.fs.tb_heater_to_cryst.Constraint() + def eq_flow_mass_comp_vapor(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + + @m.fs.tb_heater_to_cryst.Constraint() + def eq_temperature(blk): + return blk.properties_in[0].temperature == blk.properties_out[0].temperature + + @m.fs.tb_heater_to_cryst.Constraint() + def eq_pressure(blk): + return blk.properties_in[0].pressure == blk.properties_out[0].pressure + + m.fs.tb_recycle = Translator( + inlet_property_package=m.fs.properties, + outlet_property_package=m.fs.properties_nacl, + ) + + @m.fs.tb_recycle.Constraint() + def eq_flow_mass_comp(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] + ) + @m.fs.tb_recycle.Constraint() + def eq_flow_mass_comp_vapor(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + + @m.fs.tb_recycle.Constraint() + def eq_temperature(blk): + return blk.properties_in[0].temperature == blk.properties_out[0].temperature + + @m.fs.tb_recycle.Constraint() + def eq_pressure(blk): + return blk.properties_in[0].pressure == blk.properties_out[0].pressure + + + + # connections + m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.mixer.feed) + m.fs.s13 = Arc( + source=m.fs.separator.recycle, destination=m.fs.mixer.recycle) + m.fs.s11 = Arc(source=m.fs.mixer.outlet, destination=m.fs.pump.inlet) + m.fs.s02 = Arc(source=m.fs.pump.outlet, destination=m.fs.heater.cold_side_inlet) + m.fs.s03 = Arc(source=m.fs.heater.cold_side_outlet, destination=m.fs.tb_heater_to_cryst.inlet) + m.fs.s04 = Arc(source=m.fs.tb_heater_to_cryst.outlet, destination=m.fs.crystallizer.inlet) + m.fs.s05 = Arc(source=m.fs.crystallizer.outlet, destination=m.fs.tb_recycle.inlet) + m.fs.s06 = Arc(source=m.fs.tb_recycle.outlet, destination=m.fs.separator.inlet) + m.fs.s07 = Arc(source=m.fs.crystallizer.vapor, destination=m.fs.tb_vapor.inlet) + m.fs.s08 = Arc(source=m.fs.tb_vapor.outlet, destination=m.fs.compressor.inlet) + m.fs.s09 = Arc(source=m.fs.compressor.outlet, destination=m.fs.heater.hot_side_inlet) + m.fs.s10 = Arc(source=m.fs.heater.hot_side_outlet, destination=m.fs.distillate.inlet) + + + TransformationFactory("network.expand_arcs").apply_to(m) + add_costs(m) + + + + # set default property values + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e2, index=("Liq", "NaCl") + ) + m.fs.properties_vapor.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + m.fs.properties_vapor.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Vap", "H2O") + ) + + iscale.set_scaling_factor(m.fs.heater.hot.heat, 1e-3) + iscale.set_scaling_factor(m.fs.heater.cold.heat, 1e-3) + iscale.set_scaling_factor(m.fs.heater.overall_heat_transfer_coefficient, 1e-3) + iscale.set_scaling_factor(m.fs.heater.area, 1e-1) + iscale.set_scaling_factor(m.fs.compressor.control_volume.work, 1e-6) + + iscale.calculate_scaling_factors(m) + + return m + +def add_costs(m): + + m.fs.costing = WaterTAPCosting() + m.fs.crystallizer.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method_arguments={"cost_type": CrystallizerCostType.mass_basis}, + ) + m.fs.heater.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) + m.fs.mixer.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) + + m.fs.compressor.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing + ) + m.fs.pump.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method=cost_pump, + costing_method_arguments={"pump_type": "low_pressure"}, + ) + + + m.fs.costing.cost_process() + m.fs.costing.add_annual_water_production(m.fs.distillate.properties[0].flow_vol) + m.fs.costing.add_specific_energy_consumption(m.fs.distillate.properties[0].flow_vol) + m.fs.costing.add_LCOW(m.fs.distillate.properties[0].flow_vol) + +def set_operating_conditions(m): + m.fs.feed.flow_mass_phase_comp[0, "Liq", "NaCl"].fix(10.5119 ) + m.fs.feed.flow_mass_phase_comp[0, "Liq", "H2O"].fix(38.9326 ) + m.fs.crystallizer.inlet.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(1e-6) + m.fs.crystallizer.inlet.flow_mass_phase_comp[0, "Vap", "H2O"].fix(1e-6) + m.fs.feed.pressure[0].fix(101325) + m.fs.feed.temperature[0].fix(273.15 + 20) + + + m.fs.crystallizer.temperature_operating.set_value(273.15 + 50 ) + + m.fs.heater.overall_heat_transfer_coefficient.fix(2e3) + + # Fix + m.fs.crystallizer.crystal_growth_rate.fix() + m.fs.crystallizer.souders_brown_constant.fix() + m.fs.crystallizer.crystal_median_length.fix() + + + + m.fs.crystallizer.inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].set_value(13.5119 *10 ) + m.fs.crystallizer.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].set_value(38.9326 * 10 ) + m.fs.crystallizer.inlet.temperature[0].fix(273.15 + 55) + m.fs.heater.cold_side_outlet.temperature[0].set_value(273.15 + 55) + m.fs.heater.cold_side_inlet.temperature[0].set_value(273.15 + 50) + m.fs.crystallizer.inlet.pressure[0].set_value(101325) + + m.fs.compressor.pressure_ratio.fix(1.5) + m.fs.compressor.efficiency.fix(0.8) + m.fs.compressor.pressure_ratio.setub(3) + + m.fs.pump.deltaP.fix(100000) + m.fs.pump.efficiency_pump.fix(0.8) + + m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(2) + + + + print("DOF:", degrees_of_freedom(m.fs)) + +def solve(blk, solver=None, tee=True): + if solver is None: + solver = get_solver() + results = solver.solve(blk, tee=tee) + if not check_optimal_termination(results): + results = solver.solve(blk, tee=tee) + return results + + +def initialize_system(m, solver=None, verbose=True): + if solver is None: + solver = get_solver() + + # initialize feed block + m.fs.feed.initialize() + m.fs.crystallizer.initialize() + + propagate_state(m.fs.s07) + propagate_state(m.fs.s08) + m.fs.compressor.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.compressor.inlet.flow_mass_phase_comp[0,"Vap", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"].value + m.fs.compressor.inlet.temperature[0] = m.fs.crystallizer.vapor.temperature[0].value + m.fs.compressor.inlet.pressure[0] = m.fs.crystallizer.vapor.pressure[0].value + m.fs.compressor.initialize() + + propagate_state(m.fs.s01) + + propagate_state(m.fs.s05) + propagate_state(m.fs.s06) + + m.fs.separator.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.separator.inlet.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value + m.fs.separator.inlet.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value + m.fs.separator.inlet.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + m.fs.separator.split_fraction[0, "purge"].fix(0.1) + + m.fs.separator.initialize() + m.fs.separator.split_fraction[0, "purge"].unfix() + + + propagate_state(m.fs.s13) + m.fs.mixer.initialize() + m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + print(f"DOF: {degrees_of_freedom(m)}") + + propagate_state(m.fs.s11) + m.fs.pump.initialize() + + propagate_state(m.fs.s02) + propagate_state(m.fs.s09) + m.fs.heater.initialize() + m.fs.heater.cold_side_inlet.unfix() + m.fs.heater.hot_side_inlet.unfix() + print(f"DOF heater: {degrees_of_freedom(m)}") + + propagate_state(m.fs.s03) + propagate_state(m.fs.s04) + m.fs.crystallizer.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.heater.cold_side_outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.crystallizer.inlet.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.heater.cold_side_outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value + m.fs.crystallizer.inlet.temperature[0] = m.fs.heater.cold_side_outlet.temperature[0].value + m.fs.crystallizer.inlet.pressure[0] = m.fs.heater.cold_side_outlet.pressure[0].value + m.fs.crystallizer.initialize() + print(f"DOF: {degrees_of_freedom(m)}") + propagate_state(m.fs.s10) + m.fs.distillate.initialize() + print("DOF final:", degrees_of_freedom(m.fs)) + propagate_state(m.fs.s07) + propagate_state(m.fs.s08) + m.fs.compressor.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.compressor.inlet.flow_mass_phase_comp[0,"Vap", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"].value + m.fs.compressor.inlet.temperature[0] = m.fs.crystallizer.vapor.temperature[0].value + m.fs.compressor.inlet.pressure[0] = m.fs.crystallizer.vapor.pressure[0].value + m.fs.compressor.initialize() + m.fs.costing.initialize() + + + +def optimize_set_up(m): + m.fs.objective = Objective(expr=m.fs.costing.LCOW) + m.fs.compressor.pressure_ratio.unfix() + m.fs.compressor.pressure_ratio.setlb(1) + m.fs.compressor.pressure_ratio.setub(3) + #additional constraint + m.fs.eq_heater_temperature_rise = Constraint( + expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] <= 4) + + m.fs.crystallizer.inlet.temperature.unfix() + m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) + m.fs.crystallizer.inlet.temperature.setub(273.15 +110) + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].lb = 100 + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].set_value(500) + + +def optimize(m, solver=None): + # --solve--- + return solve(m, solver=solver) + +def display_system(m): + print("operating temperature:", m.fs.crystallizer.temperature_operating.value) + print("Levelized cost of water: %.2f $/m3" % value(m.fs.costing.LCOW)) + print("Inlet temperature heater:", m.fs.heater.cold_side_inlet.temperature[0].value) + print("compressor pressure ratio:", m.fs.compressor.pressure_ratio.value) + print("separator split recycle:", m.fs.separator.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].value) + print("separator split purge:", m.fs.separator.purge.flow_mass_phase_comp[0,"Liq", "H2O"].value) + print("distillate:", m.fs.distillate.flow_mass_phase_comp[0,"Liq", "H2O"].value) +if __name__ == "__main__": + m = main() + diff --git a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py new file mode 100644 index 0000000000..b79e256695 --- /dev/null +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -0,0 +1,508 @@ +################################################################################# +# 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 pyomo.environ import ( + ConcreteModel, + TerminationCondition, +) +from pyomo.environ import ( + ConcreteModel, + value, + Constraint, + Objective, + Var, + NonNegativeReals, + TransformationFactory, + units as pyunits, + check_optimal_termination, +) +from pyomo.network import Arc, SequentialDecomposition +from idaes.core import FlowsheetBlock +from watertap.core.solvers import get_solver +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.initialization import ( + propagate_state, +) +from pyomo.util.check_units import assert_units_consistent + +from idaes.core import FlowsheetBlock +from idaes.core.util.model_statistics import degrees_of_freedom + +import idaes.core.util.scaling as iscale +import idaes.logger as idaeslog +from watertap.core.solvers import get_solver +from idaes.core import UnitModelCostingBlock + +from watertap.property_models.unit_specific import cryst_prop_pack as props +from watertap.unit_models.crystallizer import Crystallization +from watertap.costing import WaterTAPCosting, CrystallizerCostType +from watertap.costing.unit_models.heat_exchanger import ( + cost_heat_exchanger, +) + +from io import StringIO +from pyomo.util.infeasible import ( + log_infeasible_constraints, +) +from watertap.unit_models.steam_heater_0D import SteamHeater0D, Mode +from idaes.models.unit_models import Heater, Separator, Mixer, Product, Feed +from idaes.models.unit_models.separator import SplittingType +from watertap.unit_models.pressure_changer import Pump +from watertap.costing.unit_models.pump import ( + cost_pump, +) +from idaes.models.unit_models.mixer import MomentumMixingType +from watertap.core.solvers import get_solver +from watertap.unit_models.mvc.components.lmtd_chen_callback import ( + delta_temperature_chen_callback, +) +from idaes.models.unit_models.heat_exchanger import ( + HeatExchangerFlowPattern, +) +from watertap.costing.unit_models.heater_chiller import ( + cost_heater_chiller, +) +from watertap.property_models.unit_specific import cryst_prop_pack as props +import watertap.property_models.water_prop_pack as props_w +import watertap.property_models.NaCl_T_dep_prop_pack as props_nacl +from pyomo.common.log import LoggingIntercept +from idaes.models.unit_models.translator import Translator +import logging +from watertap.core.util.model_debug_mode import activate +activate() + + + +__author__ = "Elmira Shamlou" + + +def main(): + solver = get_solver() + m = build() + set_operating_conditions(m) + initialize_system(m, solver=solver) + + optimize_set_up(m) + + #interval_initializer(m) + solve(m, solver=solver) + + #print("\n***---optimization results---***") + display_system(m) + #display_design(m) + #display_state(m) + + return m + + +def build(): + + # flowsheet set up + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + # properties + m.fs.properties_vapor = props_w.WaterParameterBlock() + m.fs.properties = props.NaClParameterBlock() + m.fs.properties_nacl = props_nacl.NaClParameterBlock() + + # Control volume flow blocks + m.fs.feed = Feed(property_package= m.fs.properties_nacl) + m.fs.distillate = Product(property_package= m.fs.properties_vapor) + + + + # overall process water recovery (evaporation rate/solid content) + + + # unit models: steam heater, mixer, pump, crystalizer + + m.fs.crystallizer = Crystallization(property_package=m.fs.properties) + m.fs.heater = SteamHeater0D( + hot_side_name="hot", + cold_side_name="cold", + hot={ + "property_package": m.fs.properties_vapor, + }, + cold={ + "property_package": m.fs.properties_nacl, + }, + delta_temperature_callback=delta_temperature_chen_callback, + flow_pattern=HeatExchangerFlowPattern.countercurrent, + mode=Mode.HEATER, + ) + + m.fs.condenser = SteamHeater0D( + hot_side_name="hot", + cold_side_name="cold", + hot={ + "property_package": m.fs.properties_vapor, + }, + cold={ + "property_package": m.fs.properties_nacl, + }, + delta_temperature_callback=delta_temperature_chen_callback, + flow_pattern=HeatExchangerFlowPattern.countercurrent, + mode=Mode.CONDENSER, + estimate_cooling_water=True + ) + + m.fs.chiller = Heater( + property_package=m.fs.properties_nacl, has_pressure_change=False + ) + + m.fs.mixer = Mixer( + property_package=m.fs.properties_nacl, + momentum_mixing_type=MomentumMixingType.equality, + inlet_list=["feed", "recycle"], + ) + m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + + m.fs.pump = Pump(property_package=m.fs.properties_nacl) + + m.fs.separator = Separator( + property_package=m.fs.properties_nacl, + outlet_list=["recycle", "purge"], + split_basis=SplittingType.totalFlow, + ) + + m.fs.tb_recycle = Translator( + inlet_property_package=m.fs.properties, + outlet_property_package=m.fs.properties_nacl, + ) + + m.fs.tb_vapor = Translator( + inlet_property_package=m.fs.properties, + outlet_property_package=m.fs.properties_vapor + ) + + @m.fs.tb_vapor.Constraint() + def eq_flow_mass_comp(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Vap", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Vap", "H2O"] + ) + @m.fs.tb_vapor.Constraint() + def eq_flow_mass_comp_vapor(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] + ) + + @m.fs.tb_vapor.Constraint() + def eq_temperature(blk): + return blk.properties_in[0].temperature == blk.properties_out[0].temperature + + @m.fs.tb_vapor.Constraint() + def eq_pressure(blk): + return blk.properties_in[0].pressure == blk.properties_out[0].pressure + + @m.fs.tb_recycle.Constraint() + def eq_flow_mass_comp(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] + ) + @m.fs.tb_recycle.Constraint() + def eq_flow_mass_comp_solute(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + + @m.fs.tb_recycle.Constraint() + def eq_temperature(blk): + return blk.properties_in[0].temperature == blk.properties_out[0].temperature + + @m.fs.tb_recycle.Constraint() + def eq_pressure(blk): + return blk.properties_in[0].pressure == blk.properties_out[0].pressure + + m.fs.tb_heater_to_cryst = Translator( + inlet_property_package=m.fs.properties_nacl, + outlet_property_package=m.fs.properties, + ) + + @m.fs.tb_heater_to_cryst.Constraint() + def eq_flow_mass_comp(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] + ) + @m.fs.tb_heater_to_cryst.Constraint() + def eq_flow_mass_comp_solute(blk): + return ( + blk.properties_in[0].flow_mass_phase_comp["Liq", "NaCl"] + == blk.properties_out[0].flow_mass_phase_comp["Liq", "NaCl"] + ) + + @m.fs.tb_heater_to_cryst.Constraint() + def eq_temperature(blk): + return blk.properties_in[0].temperature == blk.properties_out[0].temperature + + @m.fs.tb_heater_to_cryst.Constraint() + def eq_pressure(blk): + return blk.properties_in[0].pressure == blk.properties_out[0].pressure + + + m.fs.eq_heater_temperature_rise = Constraint( + expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] == 4 * pyunits.K) + m.fs.eq_chiller = Constraint( + expr=m.fs.chiller.control_volume.properties_out[0].temperature == m.fs.condenser.cold_side_inlet.temperature[0]) + + # performance expressions + + # connections + m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.mixer.feed) + m.fs.s02 = Arc( + source=m.fs.separator.recycle, destination=m.fs.mixer.recycle) + m.fs.s03 = Arc(source=m.fs.mixer.outlet, destination=m.fs.pump.inlet) + m.fs.s04 = Arc(source=m.fs.pump.outlet, destination=m.fs.heater.cold_side_inlet) + m.fs.s05 = Arc(source=m.fs.heater.cold_side_outlet, destination=m.fs.tb_heater_to_cryst.inlet) + m.fs.s06 = Arc(source=m.fs.tb_heater_to_cryst.outlet, destination=m.fs.crystallizer.inlet) + + m.fs.s07 = Arc(source=m.fs.crystallizer.outlet, destination=m.fs.tb_recycle.inlet) + m.fs.s08 = Arc(source=m.fs.tb_recycle.outlet, destination=m.fs.separator.inlet) + + m.fs.s09 = Arc(source=m.fs.crystallizer.vapor, destination=m.fs.tb_vapor.inlet) + m.fs.s10 = Arc(source=m.fs.tb_vapor.outlet, destination=m.fs.condenser.hot_side_inlet) + m.fs.s11 = Arc(source=m.fs.condenser.hot_side_outlet, destination=m.fs.distillate.inlet) + m.fs.s12 = Arc(source=m.fs.condenser.cold_side_outlet, destination=m.fs.chiller.inlet) + + + + + + TransformationFactory("network.expand_arcs").apply_to(m) + add_costs(m) + + + + # set default property values + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mass_phase_comp", 1e2, index=("Liq", "NaCl") + ) + m.fs.properties_vapor.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Liq", "H2O") + ) + m.fs.properties_vapor.set_default_scaling( + "flow_mass_phase_comp", 1, index=("Vap", "H2O") + ) + + iscale.set_scaling_factor(m.fs.heater.hot.heat, 1e-3) + iscale.set_scaling_factor(m.fs.heater.cold.heat, 1e-3) + iscale.set_scaling_factor(m.fs.heater.overall_heat_transfer_coefficient, 1e-3) + iscale.set_scaling_factor(m.fs.heater.area, 1e-1) + + iscale.set_scaling_factor(m.fs.condenser.hot.heat, 1e-3) + iscale.set_scaling_factor(m.fs.condenser.cold.heat, 1e-3) + iscale.set_scaling_factor(m.fs.condenser.overall_heat_transfer_coefficient, 1e-3) + iscale.set_scaling_factor(m.fs.condenser.area, 1e-1) + iscale.set_scaling_factor(m.fs.chiller.control_volume.heat, 1e-3) + + iscale.calculate_scaling_factors(m) + + return m + +def add_costs(m): + + m.fs.costing = WaterTAPCosting() + m.fs.crystallizer.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method_arguments={"cost_type": CrystallizerCostType.mass_basis}, +) + m.fs.heater.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing, + costing_method=cost_heat_exchanger, + costing_method_arguments={"cost_steam_flow": True},) + + m.fs.condenser.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing, + costing_method=cost_heat_exchanger) + + m.fs.chiller.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method=cost_heater_chiller, + costing_method_arguments={"HC_type": "chiller"}, + ) + + m.fs.mixer.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) + m.fs.pump.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method=cost_pump, + costing_method_arguments={"pump_type": "low_pressure"}, + ) + + m.fs.costing.cost_process() + m.fs.costing.add_annual_water_production(m.fs.distillate.properties[0].flow_vol) + m.fs.costing.add_specific_energy_consumption(m.fs.distillate.properties[0].flow_vol) + m.fs.costing.add_LCOW(m.fs.distillate.properties[0].flow_vol) + +def set_operating_conditions(m): + m.fs.feed.flow_mass_phase_comp[0, "Liq", "NaCl"].fix(9.5119 ) + m.fs.feed.flow_mass_phase_comp[0, "Liq", "H2O"].fix(38.9326 ) + m.fs.tb_heater_to_cryst.outlet.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(1e-6) + m.fs.tb_heater_to_cryst.outlet.flow_mass_phase_comp[0, "Vap", "H2O"].fix(1e-6) + m.fs.feed.pressure[0].fix(101325) + m.fs.feed.temperature[0].fix(273.15 + 20) + + + m.fs.crystallizer.inlet.temperature[0].set_value(273.15 + 54 ) + #m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(5) + m.fs.crystallizer.temperature_operating.set_value(273.15 + 50 ) + m.fs.heater.overall_heat_transfer_coefficient.fix(2e3) + m.fs.heater.area.set_value(500) + + # Fix + m.fs.crystallizer.crystal_growth_rate.fix() + m.fs.crystallizer.souders_brown_constant.fix() + m.fs.crystallizer.crystal_median_length.fix() + + m.fs.heater.hot_side_inlet.flow_mass_phase_comp[0, "Vap", "H2O"].set_value(30) + m.fs.heater.hot_side_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(0) + m.fs.heater.hot_side_inlet.temperature.fix(273.15 + 170) + m.fs.heater.cold_side_outlet.temperature.fix(273.15 + 54) + + m.fs.heater.hot_side_inlet.pressure[0].fix(201325) + + m.fs.crystallizer.inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].set_value(13.5119 *10 ) + m.fs.crystallizer.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].set_value(38.9326 * 10 ) + #m.fs.crystallizer.inlet.temperature[0].set_value(273.15 + 80) + m.fs.crystallizer.inlet.pressure[0].set_value(101325) + + #m.fs.heater.area.setlb(10) + + m.fs.pump.deltaP.fix(100000) + m.fs.pump.efficiency_pump.fix(0.8) + + m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(2) + + m.fs.condenser.cold_side_inlet.pressure[0].fix(101325) + m.fs.condenser.cold_side_inlet.temperature[0].fix(273.15 + 20) + m.fs.condenser.cold_side_outlet.temperature[0].fix(273.15 + 35) + m.fs.condenser.overall_heat_transfer_coefficient.fix(2e3) + m.fs.condenser.cold_side_inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].fix( + 1e-8 + ) + m.fs.condenser.cold_side_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].set_value( + 12 + ) + m.fs.condenser.area.set_value(10) + + print("DOF final:", degrees_of_freedom(m.fs)) + +def solve(blk, solver=None, tee=True): + if solver is None: + solver = get_solver() + results = solver.solve(blk, tee=tee) + if not check_optimal_termination(results): + results = solver.solve(blk, tee=tee) + return results + + +def initialize_system(m, solver=None, verbose=True): + if solver is None: + solver = get_solver() + + # initialize feed block + m.fs.feed.initialize() + m.fs.crystallizer.initialize() + + + #separator + propagate_state(m.fs.s07) + propagate_state(m.fs.s08) + + m.fs.separator.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.separator.inlet.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value + m.fs.separator.inlet.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value + m.fs.separator.inlet.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + m.fs.separator.split_fraction[0, "purge"].fix(0.1) + + m.fs.separator.initialize() + m.fs.separator.split_fraction[0, "purge"].unfix() + + #mixer + propagate_state(m.fs.s01) + propagate_state(m.fs.s02) + + m.fs.mixer.initialize() + m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + + propagate_state(m.fs.s03) + m.fs.pump.initialize() + + propagate_state(m.fs.s04) + + m.fs.heater.initialize() + m.fs.heater.cold_side_inlet.unfix() + + propagate_state(m.fs.s05) + propagate_state(m.fs.s06) + m.fs.crystallizer.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.heater.cold_side_outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.crystallizer.inlet.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.heater.cold_side_outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value + m.fs.crystallizer.inlet.temperature[0] = m.fs.heater.cold_side_outlet.temperature[0].value + m.fs.crystallizer.inlet.pressure[0] = m.fs.heater.cold_side_outlet.pressure[0].value + + m.fs.crystallizer.initialize() + + propagate_state(m.fs.s09) + propagate_state(m.fs.s10) + m.fs.condenser.hot_side_inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.condenser.hot_side_inlet.flow_mass_phase_comp[0,"Vap", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"].value + m.fs.condenser.hot_side_inlet.temperature[0] = m.fs.crystallizer.vapor.temperature[0].value + m.fs.condenser.hot_side_inlet.pressure[0] = m.fs.crystallizer.vapor.pressure[0].value + m.fs.condenser.initialize() + m.fs.condenser.hot_side_inlet.unfix() + + propagate_state(m.fs.s11) + m.fs.distillate.initialize() + + propagate_state(m.fs.s12) + m.fs.chiller.initialize() + + + m.fs.costing.initialize() + print(f"DOF: {degrees_of_freedom(m)}") + +def optimize_set_up(m): + # add objective + m.fs.objective = Objective(expr=m.fs.costing.LCOW) + + + m.fs.heater.cold_side_outlet.temperature.unfix() + m.fs.heater.cold_side_outlet.temperature.setlb(273.15 +50) + m.fs.heater.cold_side_outlet.temperature.setub(273.15 +110) + + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].lb = 300 + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].set_value(400) + + + print("DOF final 2:", degrees_of_freedom(m.fs)) + + +def optimize(m, solver=None): + # --solve--- + return solve(m, solver=solver) +def display_system(m): + print("Crystalizer Inlet temperature:", m.fs.crystallizer.inlet.temperature[0].value) + print("operating temperature:", m.fs.crystallizer.temperature_operating.value) + print("recycle :", m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].value) + print("Levelized cost of water: %.2f $/m3" % value(m.fs.costing.LCOW)) + print("steam:", m.fs.heater.hot_side_inlet.flow_mass_phase_comp[0,"Vap", "H2O"].value) + print("heater area:", m.fs.heater.area.value) + + + + +if __name__ == "__main__": + m = main() diff --git a/watertap/unit_models/Crystallizer_revised.py b/watertap/unit_models/Crystallizer_revised.py new file mode 100644 index 0000000000..580a3bbe33 --- /dev/null +++ b/watertap/unit_models/Crystallizer_revised.py @@ -0,0 +1,827 @@ +################################################################################# +# 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, + check_optimal_termination, + Param, + Constraint, + Suffix, + units as pyunits, +) +from pyomo.common.config import ConfigBlock, ConfigValue, In + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + UnitModelBlockData, + useDefault, +) +from watertap.core.solvers import get_solver +from idaes.core.util.tables import create_stream_table_dataframe +from idaes.core.util.constants import Constants +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 import InitializationMixin +from watertap.costing.unit_models.crystallizer import cost_crystallizer +from watertap.core.util.initialization import interval_initializer + +_log = idaeslog.getLogger(__name__) + +__author__ = "Oluwamayowa Amusat" + + +# when using this file the name "Filtration" is what is imported +@declare_process_block_class("Crystallization") +class CrystallizationData(InitializationMixin, UnitModelBlockData): + """ + Zero order crystallization model + """ + + # CONFIG are options for the unit model, this simple model only has the mandatory config options + 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_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.}""", + ), + ) + + def build(self): + super().build() + + solvent_set = self.config.property_package.solvent_set + solute_set = self.config.property_package.solute_set + + # this creates blank scaling factors, which are populated later + self.scaling_factor = Suffix(direction=Suffix.EXPORT) + + # Next, get the base units of measurement from the property definition + units_meta = self.config.property_package.get_metadata().get_derived_units + + # Add unit variables + """ + self.approach_temperature_heat_exchanger = Param( + initialize=4, + units=pyunits.K, + doc="Maximum temperature difference between inlet and outlet of a crystallizer heat exchanger.\ + Lewis et al. suggests 1-2 degC but use 5degC in example; Tavare example used 4 degC.\ + Default is 4 degC", + ) + """ + + # ====== Crystallizer sizing parameters ================= # + self.dimensionless_crystal_length = Param( + initialize=3.67, # Parameter from population balance modeling for median crystal size + units=pyunits.dimensionless, + ) + + self.crystal_median_length = Var( + initialize=0.5e-3, # From Mersmann et al., Tavare et al. example + bounds=( + 0.2e-3, + 0.6e-3, + ), # Limits for FC crystallizers based on Bermingham et al. + units=pyunits.m, + doc="Desired median crystal size, m", + ) + + self.crystal_growth_rate = Var( + initialize=3.7e-8, # From Mersmann et al. for NaCl. Perry has values between 0.5e-8 to 13e-8 for NaCl + bounds=(1e-9, 1e-6), # Based on Mersmann and Kind diagram. + units=pyunits.m / pyunits.s, + doc="Crystal growth rate, m/s", + ) + + self.souders_brown_constant = Var( + initialize=0.04, + units=pyunits.m / pyunits.s, + doc="Constant for Souders-Brown equation, set at 0.04 m/s based on Dutta et al. \ + Lewis et al suggests 0.024 m/s, while Tavare suggests about 0.06 m/s ", + ) + + # ====== Model variables ================= # + self.crystallization_yield = Var( + solute_set, + initialize=0.5, + bounds=(0.0, 1), + units=pyunits.dimensionless, + doc="Crystallizer solids yield", + ) + + self.product_volumetric_solids_fraction = Var( + initialize=0.25, + bounds=(0.0, 1), + units=pyunits.dimensionless, + doc="Volumetric fraction of solids in slurry product (i.e. solids-liquid mixture).", + ) + + self.temperature_operating = Var( + initialize=298.15, + bounds=(273, 1000), + units=pyunits.K, + doc="Crystallizer operating temperature: boiling point of the solution.", + ) + + self.pressure_operating = Var( + initialize=1e3, + bounds=(0.001, 2e6), + units=pyunits.Pa, + doc="Operating pressure of the crystallizer.", + ) + + self.dens_mass_magma = Var( + initialize=250, + bounds=(0.01, 5000), + units=pyunits.kg / pyunits.m**3, + doc="Magma density, i.e. mass of crystals per unit volume of suspension", + ) + + self.dens_mass_slurry = Var( + initialize=1000, + bounds=(1, 5000), + units=pyunits.kg / pyunits.m**3, + doc="Suspension density, i.e. density of solid-liquid mixture before separation", + ) + """ + self.work_mechanical = Var( + self.flowsheet().config.time, + initialize=1e5, + bounds=(-5e6, 5e6), + units=pyunits.kJ / pyunits.s, + doc="Crystallizer thermal energy requirement", + ) + """ + self.diameter_crystallizer = Var( + initialize=3, + bounds=(0, 25), + units=pyunits.m, + doc="Diameter of crystallizer", + ) + + self.height_slurry = Var( + initialize=3, + bounds=(0, 25), + units=pyunits.m, + doc="Slurry height in crystallizer", + ) + + self.height_crystallizer = Var( + initialize=3, bounds=(0, 30), units=pyunits.m, doc="Crystallizer height" + ) + + + """ + self.magma_circulation_flow_vol = Var( + initialize=1, + bounds=(0, 100), + units=pyunits.m**3 / pyunits.s, + doc="Minimum circulation flow rate through crystallizer heat exchanger", + ) + """ + + self.relative_supersaturation = Var( + solute_set, initialize=0.1, bounds=(0, 100), units=pyunits.dimensionless + ) + + self.t_res = Var( + initialize=1, + bounds=(0, 10), + units=pyunits.hr, + doc="Residence time in crystallizer", + ) + + self.volume_suspension = Var( + initialize=1, + bounds=(0, None), + units=pyunits.m**3, + doc="Crystallizer minimum active volume, i.e. volume of liquid-solid suspension", + ) + + # Add state blocks for inlet, outlet, and waste + # These include the state variables and any other properties on demand + # Add inlet block + 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"] = True # inlet block is an inlet + self.properties_in = self.config.property_package.state_block_class( + self.flowsheet().config.time, doc="Material properties of inlet", **tmp_dict + ) + + # Add outlet and waste block + tmp_dict["defined_state"] = False # outlet and waste block is not an inlet + self.properties_out = self.config.property_package.state_block_class( + self.flowsheet().config.time, + doc="Material properties of liquid outlet", + **tmp_dict, + ) + + self.properties_solids = self.config.property_package.state_block_class( + self.flowsheet().config.time, + doc="Material properties of solid crystals at outlet", + **tmp_dict, + ) + + self.properties_vapor = self.config.property_package.state_block_class( + self.flowsheet().config.time, + doc="Material properties of water vapour at outlet", + **tmp_dict, + ) + + # Add ports - oftentimes users interact with these rather than the state blocks + self.add_port(name="inlet", block=self.properties_in) + self.add_port(name="outlet", block=self.properties_out) + self.add_port(name="solids", block=self.properties_solids) + self.add_port(name="vapor", block=self.properties_vapor) + + # Add constraints + # 1. Material balances + @self.Constraint( + self.config.property_package.component_list, + doc="Mass balance for components", + ) + def eq_mass_balance_constraints(b, j): + return sum( + b.properties_in[0].flow_mass_phase_comp[p, j] + for p in self.config.property_package.phase_list + if (p, j) in b.properties_in[0].phase_component_set + ) == sum( + b.properties_out[0].flow_mass_phase_comp[p, j] + for p in self.config.property_package.phase_list + if (p, j) in b.properties_out[0].phase_component_set + ) + sum( + b.properties_vapor[0].flow_mass_phase_comp[p, j] + for p in self.config.property_package.phase_list + if (p, j) in b.properties_vapor[0].phase_component_set + ) + sum( + b.properties_solids[0].flow_mass_phase_comp[p, j] + for p in self.config.property_package.phase_list + if (p, j) in b.properties_solids[0].phase_component_set + ) + + # 2. Constraint on outlet liquid composition based on solubility requirements + @self.Constraint( + self.config.property_package.component_list, + doc="Solubility vs mass fraction constraint", + ) + def eq_solubility_massfrac_equality_constraint(b, j): + if j in solute_set: + return ( + b.properties_out[0].mass_frac_phase_comp["Liq", j] + - b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] + == 0 + ) + else: + return Constraint.Skip + + """ + # 3. Performance equations + # (a) based on yield + @self.Constraint( + self.config.property_package.component_list, + doc="Component salt yield equation", + ) + def eq_removal_balance(b, j): + if j in solvent_set: + return Constraint.Skip + else: + return ( + b.properties_in[0].flow_mass_phase_comp["Liq", j] + * b.crystallization_yield[j] + == b.properties_in[0].flow_mass_phase_comp["Liq", j] + - b.properties_out[0].flow_mass_phase_comp["Liq", j] + ) + + # (b) Volumetric fraction constraint + @self.Constraint(doc="Solid volumetric fraction in outlet: constraint, 1-E") + def eq_vol_fraction_solids(b): + return self.product_volumetric_solids_fraction * ( + b.properties_solids[0].flow_vol + b.properties_out[0].flow_vol + )== b.properties_solids[ + 0 + ].flow_vol + + # (c) Magma density constraint + @self.Constraint(doc="Slurry magma density") + def eq_dens_magma(b): + return ( + self.dens_mass_magma + == b.properties_solids[0].dens_mass_solute["Sol"] + * self.product_volumetric_solids_fraction + ) + + # (e) Relative supersaturation + @self.Constraint( + solute_set, + doc="Relative supersaturation created via evaporation, g/g (solution)", + ) + + def eq_relative_supersaturation(b, j): + # Calculate the mass of solute j after evaporation + #### change the naming + mass_after_evap = b.properties_in[0].flow_mass_phase_comp["Liq", j] + + # Calculate the total mass of liquid phase components after evaporation + total_mass_after_evap = ( + sum( + b.properties_in[0].flow_mass_phase_comp["Liq", k] + for k in b.config.property_package.solute_set + ) + + b.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + - b.properties_vapor[0].flow_mass_phase_comp["Vap", "H2O"] + ) + + return ( + b.relative_supersaturation[j] * b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] * total_mass_after_evap + == ( + mass_after_evap + - b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] * total_mass_after_evap + )) +""" + # (d) Operating pressure constraint + @self.Constraint(doc="Operating pressure constraint") + def eq_operating_pressure_constraint(b): + return self.pressure_operating == b.properties_out[0].pressure_sat + + + + + # 4. Fix flows of empty solid, liquid and vapour streams + # (i) Fix solids: liquid and vapour flows must be zero + for p, j in self.properties_solids[0].phase_component_set: + if p != "Sol": + self.properties_solids[0].flow_mass_phase_comp[p, j].fix(1e-8) + + # (ii) Fix liquids: solid and vapour flows must be zero + for p, j in self.properties_out[0].phase_component_set: + if p != "Liq": + self.properties_out[0].flow_mass_phase_comp[p, j].fix(1e-8) + + # (iii) Fix vapor: solid and vapour flows must be zero. + for p, j in self.properties_vapor[0].phase_component_set: + if p != "Vap": + self.properties_vapor[0].flow_mass_phase_comp[p, j].fix(1e-8) + + # 5. Add an energy balance for the system + ## (iii) Enthalpy balance: based on Lewis et al. Enthalpy is exothermic, and the hC in property package is -ve + @self.Constraint(doc="Enthalpy balance over crystallization system") + def eq_enthalpy_balance(b): + return ( + b.properties_in[0].enth_flow + - b.properties_out[0].enth_flow + - b.properties_vapor[0].enth_flow + - b.properties_solids[0].enth_flow + #+ self.work_mechanical[0] + - sum( + b.properties_solids[0].flow_mass_phase_comp["Sol", j] + * b.properties_solids[0].dh_crystallization_mass_comp[j] + for j in solute_set + ) + == 0 + ) + + # 6. Pressure and temperature balances - what is the pressure of the outlet solid and vapour? + # TO-DO: Figure out actual liquid and solid pressures. + @self.Constraint() + def eq_p_con1(b): + return b.properties_out[0].pressure_sat == b.properties_out[0].pressure + + @self.Constraint() + def eq_p_con2(b): + return b.properties_out[0].pressure_sat == b.properties_solids[0].pressure + + @self.Constraint() + def eq_p_con3(b): + return b.properties_vapor[0].pressure == self.pressure_operating + + @self.Constraint() + def eq_T_con1(b): + return self.temperature_operating == b.properties_solids[0].temperature + + @self.Constraint() + def eq_T_con2(b): + return self.temperature_operating == b.properties_vapor[0].temperature + + @self.Constraint() + def eq_T_con3(b): + return self.temperature_operating == b.properties_out[0].temperature + + """ + + # 7. Heat exchanger minimum circulation flow rate calculations - see Lewis et al. or Tavare et al. + @self.Constraint( + doc="Constraint on mimimum circulation rate through crystallizer heat exchanger" + ) + def eq_minimum_hex_circulation_rate_constraint(b): + dens_cp_avg = self.approach_temperature_heat_exchanger * ( + b.product_volumetric_solids_fraction + * b.properties_solids[0].dens_mass_solute["Sol"] + * b.properties_solids[0].cp_mass_solute["Sol"] + + (1 - b.product_volumetric_solids_fraction) + * b.properties_out[0].dens_mass_phase["Liq"] + * b.properties_out[0].cp_mass_phase["Liq"] + ) + return b.magma_circulation_flow_vol * dens_cp_avg == pyunits.convert( + b.work_mechanical[0], to_units=pyunits.J / pyunits.s + ) + + + # 8. Suspension density + @self.Constraint(doc="Slurry density calculation") + def eq_dens_mass_slurry(b): + return ( + self.dens_mass_slurry + == b.product_volumetric_solids_fraction + * b.properties_solids[0].dens_mass_solute["Sol"] + + (1 - b.product_volumetric_solids_fraction) + * b.properties_out[0].dens_mass_phase["Liq"] + ) + + # 9. Residence time calculation + @self.Constraint(doc="Residence time") + def eq_residence_time(b): + return b.t_res == b.crystal_median_length / ( + b.dimensionless_crystal_length + * pyunits.convert( + b.crystal_growth_rate, to_units=pyunits.m / pyunits.hr + ) + ) + + # 10. Suspension volume calculation + @self.Constraint(doc="Suspension volume") + def eq_suspension_volume(b): + return b.volume_suspension == ( + b.properties_solids[0].flow_vol + b.properties_out[0].flow_vol + ) * pyunits.convert(b.t_res, to_units=pyunits.s) + + # 11. Minimum diameter of evaporation zone + @self.Expression(doc="maximum allowable vapour linear velocity in m/s") + def eq_max_allowable_velocity(b): + return ( + b.souders_brown_constant + * ( + b.properties_out[0].dens_mass_phase["Liq"] + / b.properties_vapor[0].dens_mass_solvent["Vap"] + ) + ** 0.5 + ) + + @self.Constraint( + doc="Crystallizer diameter (based on minimum diameter of evaporation zone)" + ) + def eq_vapor_head_diameter_constraint(b): + return ( + self.diameter_crystallizer + == ( + 4 + * b.properties_vapor[0].flow_vol_phase["Vap"] + / (Constants.pi * b.eq_max_allowable_velocity) + ) + ** 0.5 + ) + + # 12. Minimum crystallizer height + @self.Constraint(doc="Slurry height based on crystallizer diameter") + def eq_slurry_height_constraint(b): + return self.height_slurry * ( + Constants.pi * b.diameter_crystallizer**2 + )== 4 * b.volume_suspension + + @self.Expression( + doc="Recommended height of vapor space (0.75*D) based on Tavares et. al." + ) + def eq_vapor_space_height(b): + return 0.75 * b.diameter_crystallizer + + @self.Expression( + doc="Height to diameter ratio constraint for evaporative crystallizers (Wilson et. al.)" + ) + def eq_minimum_height_diameter_ratio(b): + return 1.5 * b.diameter_crystallizer + + @self.Constraint(doc="Crystallizer height") + def eq_crystallizer_height_constraint(b): + # Height is max(). Manual smooth max implementation used here: max(a,b) = 0.5(a + b + |a-b|) + a = b.eq_vapor_space_height + b.height_slurry + b = b.eq_minimum_height_diameter_ratio + eps = 1e-20 * pyunits.m + return self.height_crystallizer == 0.5 * ( + a + b + ((a - b) ** 2 + eps**2) ** 0.5 + ) + + """ + + 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) + + # --------------------------------------------------------------------- + # Initialize holdup block + 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.") + # --------------------------------------------------------------------- + # Initialize other state blocks + # Set state_args from inlet state + 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 == "Liq" or p == "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, + ) + init_log.info_high("Initialization Step 2 Complete.") + + interval_initializer(self) + # --------------------------------------------------------------------- + # Solve unit + 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))) + # --------------------------------------------------------------------- + # Release Inlet state + self.properties_in.release_state(flags, outlvl=outlvl) + init_log.info("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() + + iscale.set_scaling_factor( + self.crystal_growth_rate, 1e7 + ) # growth rates typically of order 1e-7 to 1e-9 m/s + iscale.set_scaling_factor( + self.crystal_median_length, 1e3 + ) # Crystal lengths typically in mm + iscale.set_scaling_factor( + self.souders_brown_constant, 1e2 + ) # Typical values are 0.0244, 0.04 and 0.06 + iscale.set_scaling_factor( + self.diameter_crystallizer, 1 + ) # Crystallizer diameters typically up to about 20 m + iscale.set_scaling_factor( + self.height_crystallizer, 1 + ) # H/D ratio maximum is about 1.5, so same scaling as diameter + iscale.set_scaling_factor(self.height_slurry, 1) # Same scaling as diameter + #iscale.set_scaling_factor(self.magma_circulation_flow_vol, 1) + iscale.set_scaling_factor(self.relative_supersaturation, 10) + iscale.set_scaling_factor(self.t_res, 1) # Residence time is in hours + iscale.set_scaling_factor( + self.volume_suspension, 0.1 + ) # Suspension volume usually in tens to hundreds range + iscale.set_scaling_factor( + self.crystallization_yield, 1 + ) # Yield is between 0 and 1, usually in the 10-60% range + iscale.set_scaling_factor(self.product_volumetric_solids_fraction, 10) + iscale.set_scaling_factor( + self.temperature_operating, + iscale.get_scaling_factor(self.properties_in[0].temperature), + ) + iscale.set_scaling_factor(self.pressure_operating, 1e-3) + iscale.set_scaling_factor( + self.dens_mass_magma, 1e-3 + ) # scaling factor of dens_mass_phase['Liq'] + iscale.set_scaling_factor( + self.dens_mass_slurry, 1e-3 + ) # scaling factor of dens_mass_phase['Liq'] + #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"]), + #) + iscale.set_scaling_factor(self.properties_in[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_out[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_vapor[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_solids[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_in[0.0].pressure, 1e3) + + + + # transforming constraints + for ind, c in self.eq_T_con1.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].temperature) + iscale.constraint_scaling_transform(c, sf) + + for ind, c in self.eq_T_con2.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].temperature) + iscale.constraint_scaling_transform(c, sf) + + for ind, c in self.eq_T_con3.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].temperature) + iscale.constraint_scaling_transform(c, sf) + + for ind, c in self.eq_p_con1.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].pressure) + iscale.constraint_scaling_transform(c, sf) + + for ind, c in self.eq_p_con2.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].pressure) + iscale.constraint_scaling_transform(c, sf) + + for ind, c in self.eq_p_con3.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].pressure) + iscale.constraint_scaling_transform(c, sf) + + for j, c in self.eq_mass_balance_constraints.items(): + sf = iscale.get_scaling_factor( + self.properties_in[0].flow_mass_phase_comp["Liq", j] + ) + iscale.constraint_scaling_transform(c, sf) + + for j, c in self.eq_solubility_massfrac_equality_constraint.items(): + iscale.constraint_scaling_transform(c, 1e0) + + #for j, c in self.eq_dens_magma.items(): + # iscale.constraint_scaling_transform( + # c, iscale.get_scaling_factor(self.dens_mass_magma) + #) + + #for j, c in self.eq_removal_balance.items(): + # sf = iscale.get_scaling_factor( + # self.properties_in[0].flow_mass_phase_comp["Liq", j] + # ) + # iscale.constraint_scaling_transform(c, sf) + + for ind, c in self.eq_enthalpy_balance.items(): + sf = iscale.get_scaling_factor(self.properties_in[0].enth_flow) + iscale.constraint_scaling_transform(c, sf) + + + def _get_stream_table_contents(self, time_point=0): + return create_stream_table_dataframe( + { + "Feed Inlet": self.inlet, + "Liquid Outlet": self.outlet, + "Vapor Outlet": self.vapor, + "Solid Outlet": self.solids, + }, + time_point=time_point, + ) + + def _get_performance_contents(self, time_point=0): + var_dict = {} + var_dict["Operating Temperature (K)"] = self.temperature_operating + var_dict["Operating Pressure (Pa)"] = self.pressure_operating + var_dict["Magma density of solution (Kg/m**3)"] = self.dens_mass_magma + var_dict["Slurry density (Kg/m3)"] = self.dens_mass_slurry + #var_dict["Heat requirement"] = self.work_mechanical[time_point] + var_dict["Crystallizer diameter (m)"] = self.diameter_crystallizer + #var_dict["Magma circulation flow rate (m**3/s)"] = ( + # self.magma_circulation_flow_vol + # ) + var_dict["Vol. frac. of solids in suspension, 1-E"] = ( + self.product_volumetric_solids_fraction + ) + var_dict["Residence time"] = self.t_res + var_dict["Crystallizer minimum active volume (m**3)"] = self.volume_suspension + var_dict["Suspension height in crystallizer (m)"] = self.height_slurry + var_dict["Crystallizer height (m)"] = self.height_crystallizer + + for j in self.config.property_package.solute_set: + yield_mem_name = f"{j} yield (fraction)" + var_dict[yield_mem_name] = self.crystallization_yield[j] + supersat_mem_name = f"{j} relative supersaturation (mass fraction basis)" + var_dict[supersat_mem_name] = self.relative_supersaturation[j] + + return {"vars": var_dict} + + @property + def default_costing_method(self): + return cost_crystallizer diff --git a/watertap/unit_models/crystallizer.py b/watertap/unit_models/crystallizer.py index 6ed2ccaa39..580a3bbe33 100644 --- a/watertap/unit_models/crystallizer.py +++ b/watertap/unit_models/crystallizer.py @@ -40,8 +40,8 @@ import idaes.logger as idaeslog from watertap.core import InitializationMixin -from watertap.core.util.initialization import interval_initializer from watertap.costing.unit_models.crystallizer import cost_crystallizer +from watertap.core.util.initialization import interval_initializer _log = idaeslog.getLogger(__name__) @@ -122,7 +122,7 @@ def build(self): units_meta = self.config.property_package.get_metadata().get_derived_units # Add unit variables - + """ self.approach_temperature_heat_exchanger = Param( initialize=4, units=pyunits.K, @@ -130,6 +130,7 @@ def build(self): Lewis et al. suggests 1-2 degC but use 5degC in example; Tavare example used 4 degC.\ Default is 4 degC", ) + """ # ====== Crystallizer sizing parameters ================= # self.dimensionless_crystal_length = Param( @@ -186,14 +187,14 @@ def build(self): self.pressure_operating = Var( initialize=1e3, - bounds=(0.001, 1e6), + bounds=(0.001, 2e6), units=pyunits.Pa, doc="Operating pressure of the crystallizer.", ) self.dens_mass_magma = Var( initialize=250, - bounds=(1, 5000), + bounds=(0.01, 5000), units=pyunits.kg / pyunits.m**3, doc="Magma density, i.e. mass of crystals per unit volume of suspension", ) @@ -204,15 +205,15 @@ def build(self): units=pyunits.kg / pyunits.m**3, doc="Suspension density, i.e. density of solid-liquid mixture before separation", ) - - self.work_mechanical = Var( - self.flowsheet().config.time, - initialize=1e5, - bounds=(-5e6, 5e6), - units=pyunits.kJ / pyunits.s, - doc="Crystallizer thermal energy requirement", - ) - + """ + self.work_mechanical = Var( + self.flowsheet().config.time, + initialize=1e5, + bounds=(-5e6, 5e6), + units=pyunits.kJ / pyunits.s, + doc="Crystallizer thermal energy requirement", + ) + """ self.diameter_crystallizer = Var( initialize=3, bounds=(0, 25), @@ -228,15 +229,18 @@ def build(self): ) self.height_crystallizer = Var( - initialize=3, bounds=(0, 25), units=pyunits.m, doc="Crystallizer height" + initialize=3, bounds=(0, 30), units=pyunits.m, doc="Crystallizer height" ) - self.magma_circulation_flow_vol = Var( - initialize=1, - bounds=(0, 100), - units=pyunits.m**3 / pyunits.s, - doc="Minimum circulation flow rate through crystallizer heat exchanger", - ) + + """ + self.magma_circulation_flow_vol = Var( + initialize=1, + bounds=(0, 100), + units=pyunits.m**3 / pyunits.s, + doc="Minimum circulation flow rate through crystallizer heat exchanger", + ) + """ self.relative_supersaturation = Var( solute_set, initialize=0.1, bounds=(0, 100), units=pyunits.dimensionless @@ -333,31 +337,32 @@ def eq_solubility_massfrac_equality_constraint(b, j): else: return Constraint.Skip - # 3. Performance equations - # (a) based on yield - @self.Constraint( - self.config.property_package.component_list, - doc="Component salt yield equation", - ) - def eq_removal_balance(b, j): - if j in solvent_set: - return Constraint.Skip - else: - return ( - b.properties_in[0].flow_mass_phase_comp["Liq", j] - * b.crystallization_yield[j] - == b.properties_in[0].flow_mass_phase_comp["Liq", j] - - b.properties_out[0].flow_mass_phase_comp["Liq", j] - ) - + """ + # 3. Performance equations + # (a) based on yield + @self.Constraint( + self.config.property_package.component_list, + doc="Component salt yield equation", + ) + def eq_removal_balance(b, j): + if j in solvent_set: + return Constraint.Skip + else: + return ( + b.properties_in[0].flow_mass_phase_comp["Liq", j] + * b.crystallization_yield[j] + == b.properties_in[0].flow_mass_phase_comp["Liq", j] + - b.properties_out[0].flow_mass_phase_comp["Liq", j] + ) + # (b) Volumetric fraction constraint @self.Constraint(doc="Solid volumetric fraction in outlet: constraint, 1-E") def eq_vol_fraction_solids(b): - return self.product_volumetric_solids_fraction == b.properties_solids[ - 0 - ].flow_vol / ( + return self.product_volumetric_solids_fraction * ( b.properties_solids[0].flow_vol + b.properties_out[0].flow_vol - ) + )== b.properties_solids[ + 0 + ].flow_vol # (c) Magma density constraint @self.Constraint(doc="Slurry magma density") @@ -368,37 +373,41 @@ def eq_dens_magma(b): * self.product_volumetric_solids_fraction ) - # (d) Operating pressure constraint - @self.Constraint(doc="Operating pressure constraint") - def eq_operating_pressure_constraint(b): - return self.pressure_operating - b.properties_out[0].pressure_sat == 0 - # (e) Relative supersaturation @self.Constraint( solute_set, doc="Relative supersaturation created via evaporation, g/g (solution)", ) - def eq_relative_supersaturation(b, j): - # mass_frac_after_evap = SOLIDS IN + LIQUID IN - VAPOUR OUT - mass_frac_after_evap = b.properties_in[0].flow_mass_phase_comp["Liq", j] / ( - sum( - b.properties_in[0].flow_mass_phase_comp["Liq", k] - for k in solute_set - ) - + b.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] - - b.properties_vapor[0].flow_mass_phase_comp["Vap", "H2O"] - ) - # return (b.relative_supersaturation[j] * b.properties_out[0].solubility_mass_frac_phase_comp['Liq', j] == - # (mass_frac_after_evap - b.properties_out[0].solubility_mass_frac_phase_comp['Liq', j]) - # ) - return ( - b.relative_supersaturation[j] - == ( - mass_frac_after_evap - - b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] + + def eq_relative_supersaturation(b, j): + # Calculate the mass of solute j after evaporation + #### change the naming + mass_after_evap = b.properties_in[0].flow_mass_phase_comp["Liq", j] + + # Calculate the total mass of liquid phase components after evaporation + total_mass_after_evap = ( + sum( + b.properties_in[0].flow_mass_phase_comp["Liq", k] + for k in b.config.property_package.solute_set + ) + + b.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + - b.properties_vapor[0].flow_mass_phase_comp["Vap", "H2O"] ) - / b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] - ) + + return ( + b.relative_supersaturation[j] * b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] * total_mass_after_evap + == ( + mass_after_evap + - b.properties_out[0].solubility_mass_frac_phase_comp["Liq", j] * total_mass_after_evap + )) +""" + # (d) Operating pressure constraint + @self.Constraint(doc="Operating pressure constraint") + def eq_operating_pressure_constraint(b): + return self.pressure_operating == b.properties_out[0].pressure_sat + + + # 4. Fix flows of empty solid, liquid and vapour streams # (i) Fix solids: liquid and vapour flows must be zero @@ -425,11 +434,9 @@ def eq_enthalpy_balance(b): - b.properties_out[0].enth_flow - b.properties_vapor[0].enth_flow - b.properties_solids[0].enth_flow - + pyunits.convert( - self.work_mechanical[0], to_units=pyunits.J / pyunits.s - ) + #+ self.work_mechanical[0] - sum( - b.properties_solids[0].flow_mass_phase_comp["Sol", j] + b.properties_solids[0].flow_mass_phase_comp["Sol", j] * b.properties_solids[0].dh_crystallization_mass_comp[j] for j in solute_set ) @@ -440,11 +447,11 @@ def eq_enthalpy_balance(b): # TO-DO: Figure out actual liquid and solid pressures. @self.Constraint() def eq_p_con1(b): - return self.pressure_operating == b.properties_out[0].pressure + return b.properties_out[0].pressure_sat == b.properties_out[0].pressure @self.Constraint() def eq_p_con2(b): - return self.pressure_operating == b.properties_solids[0].pressure + return b.properties_out[0].pressure_sat == b.properties_solids[0].pressure @self.Constraint() def eq_p_con3(b): @@ -461,6 +468,8 @@ def eq_T_con2(b): @self.Constraint() def eq_T_con3(b): return self.temperature_operating == b.properties_out[0].temperature + + """ # 7. Heat exchanger minimum circulation flow rate calculations - see Lewis et al. or Tavare et al. @self.Constraint( @@ -478,6 +487,7 @@ def eq_minimum_hex_circulation_rate_constraint(b): return b.magma_circulation_flow_vol * dens_cp_avg == pyunits.convert( b.work_mechanical[0], to_units=pyunits.J / pyunits.s ) + # 8. Suspension density @self.Constraint(doc="Slurry density calculation") @@ -536,9 +546,9 @@ def eq_vapor_head_diameter_constraint(b): # 12. Minimum crystallizer height @self.Constraint(doc="Slurry height based on crystallizer diameter") def eq_slurry_height_constraint(b): - return self.height_slurry == 4 * b.volume_suspension / ( + return self.height_slurry * ( Constants.pi * b.diameter_crystallizer**2 - ) + )== 4 * b.volume_suspension @self.Expression( doc="Recommended height of vapor space (0.75*D) based on Tavares et. al." @@ -562,6 +572,8 @@ def eq_crystallizer_height_constraint(b): a + b + ((a - b) ** 2 + eps**2) ** 0.5 ) + """ + def initialize_build( self, state_args=None, @@ -667,7 +679,6 @@ def initialize_build( 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() @@ -687,7 +698,7 @@ def calculate_scaling_factors(self): self.height_crystallizer, 1 ) # H/D ratio maximum is about 1.5, so same scaling as diameter iscale.set_scaling_factor(self.height_slurry, 1) # Same scaling as diameter - iscale.set_scaling_factor(self.magma_circulation_flow_vol, 1) + #iscale.set_scaling_factor(self.magma_circulation_flow_vol, 1) iscale.set_scaling_factor(self.relative_supersaturation, 10) iscale.set_scaling_factor(self.t_res, 1) # Residence time is in hours iscale.set_scaling_factor( @@ -708,15 +719,20 @@ def calculate_scaling_factors(self): iscale.set_scaling_factor( self.dens_mass_slurry, 1e-3 ) # scaling factor of dens_mass_phase['Liq'] - kj_to_j = 1e3 - 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"]) - * kj_to_j, - ) + #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"]), + #) + iscale.set_scaling_factor(self.properties_in[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_out[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_vapor[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_solids[0].enth_flow, 1e-3) + iscale.set_scaling_factor(self.properties_in[0.0].pressure, 1e3) + + # transforming constraints for ind, c in self.eq_T_con1.items(): @@ -752,29 +768,21 @@ def calculate_scaling_factors(self): for j, c in self.eq_solubility_massfrac_equality_constraint.items(): iscale.constraint_scaling_transform(c, 1e0) - for j, c in self.eq_dens_magma.items(): - iscale.constraint_scaling_transform( - c, iscale.get_scaling_factor(self.dens_mass_magma) - ) + #for j, c in self.eq_dens_magma.items(): + # iscale.constraint_scaling_transform( + # c, iscale.get_scaling_factor(self.dens_mass_magma) + #) - for j, c in self.eq_removal_balance.items(): - sf = iscale.get_scaling_factor( - self.properties_in[0].flow_mass_phase_comp["Liq", j] - ) - iscale.constraint_scaling_transform(c, sf) + #for j, c in self.eq_removal_balance.items(): + # sf = iscale.get_scaling_factor( + # self.properties_in[0].flow_mass_phase_comp["Liq", j] + # ) + # iscale.constraint_scaling_transform(c, sf) for ind, c in self.eq_enthalpy_balance.items(): - sf = iscale.get_scaling_factor( - self.properties_out[0].flow_mass_phase_comp["Vap", "H2O"] - ) - sw = iscale.get_scaling_factor( - self.properties_out[0].enth_mass_solvent["Vap"] - ) - iscale.constraint_scaling_transform(c, sf * sw) + sf = iscale.get_scaling_factor(self.properties_in[0].enth_flow) + iscale.constraint_scaling_transform(c, sf) - for ind, c in self.eq_minimum_hex_circulation_rate_constraint.items(): - sf = iscale.get_scaling_factor(self.work_mechanical[0]) - iscale.constraint_scaling_transform(c, sf * 1e-3) def _get_stream_table_contents(self, time_point=0): return create_stream_table_dataframe( @@ -789,20 +797,22 @@ def _get_stream_table_contents(self, time_point=0): def _get_performance_contents(self, time_point=0): var_dict = {} - var_dict["Operating Temperature"] = self.temperature_operating - var_dict["Operating Pressure"] = self.pressure_operating - var_dict["Magma density of solution"] = self.dens_mass_magma - var_dict["Slurry density"] = self.dens_mass_slurry - var_dict["Heat requirement"] = self.work_mechanical[time_point] - var_dict["Crystallizer diameter"] = self.diameter_crystallizer - var_dict["Magma circulation flow rate"] = self.magma_circulation_flow_vol + var_dict["Operating Temperature (K)"] = self.temperature_operating + var_dict["Operating Pressure (Pa)"] = self.pressure_operating + var_dict["Magma density of solution (Kg/m**3)"] = self.dens_mass_magma + var_dict["Slurry density (Kg/m3)"] = self.dens_mass_slurry + #var_dict["Heat requirement"] = self.work_mechanical[time_point] + var_dict["Crystallizer diameter (m)"] = self.diameter_crystallizer + #var_dict["Magma circulation flow rate (m**3/s)"] = ( + # self.magma_circulation_flow_vol + # ) var_dict["Vol. frac. of solids in suspension, 1-E"] = ( self.product_volumetric_solids_fraction ) var_dict["Residence time"] = self.t_res - var_dict["Crystallizer minimum active volume"] = self.volume_suspension - var_dict["Suspension height in crystallizer"] = self.height_slurry - var_dict["Crystallizer height"] = self.height_crystallizer + var_dict["Crystallizer minimum active volume (m**3)"] = self.volume_suspension + var_dict["Suspension height in crystallizer (m)"] = self.height_slurry + var_dict["Crystallizer height (m)"] = self.height_crystallizer for j in self.config.property_package.solute_set: yield_mem_name = f"{j} yield (fraction)" diff --git a/watertap/unit_models/steam_heater_0D.py b/watertap/unit_models/steam_heater_0D.py index a6a5f9fe43..8f27a0de78 100644 --- a/watertap/unit_models/steam_heater_0D.py +++ b/watertap/unit_models/steam_heater_0D.py @@ -22,6 +22,7 @@ ) from enum import Enum, auto from pyomo.common.config import ConfigValue +from idaes.core.util.model_statistics import degrees_of_freedom _log = idaeslog.getLogger(__name__) @@ -93,7 +94,8 @@ def outlet_pressure_sat(b, t): b.hot_side.properties_out[t].pressure >= b.hot_side.properties_out[t].pressure_sat ) - + + def initialize_build(self, *args, **kwargs): """ Initialization routine for both heater and condenser modes. For condenser mode with cooling water estimation, the initialization is performed based on a specified design temperature rise on the cold side. @@ -108,8 +110,8 @@ def initialize_build(self, *args, **kwargs): self.hot_side_inlet.fix() self.cold_side_inlet.fix() self.hot_side_outlet.unfix() - self.cold_side_outlet.unfix() - self.area.fix() + #self.cold_side_outlet.unfix() + #self.area.fix() self.outlet_liquid_mass_balance.deactivate() self.outlet_pressure_sat.deactivate() @@ -139,30 +141,51 @@ def initialize_build(self, *args, **kwargs): self.config.mode == Mode.CONDENSER and not self.config.estimate_cooling_water ): + + self.hot_side_inlet.fix() + self.cold_side_inlet.fix() + #self.hot_side_outlet.unfix() + #self.cold_side_outlet.unfix() + #self.cold_side_outlet.temperature.fix() # condenser mode without cooling water estimation + self.outlet_liquid_mass_balance.deactivate() + self.outlet_pressure_sat.deactivate() super().initialize_build(*args, **kwargs) opt = get_solver(solver, optarg) + #self.cold_side_outlet.temperature.unfix() + + self.outlet_liquid_mass_balance.activate() + self.outlet_pressure_sat.activate() with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(self, tee=slc.tee) init_log.info("Initialization Complete {}".format(idaeslog.condition(res))) elif self.config.mode == Mode.CONDENSER and self.config.estimate_cooling_water: + print(f"DOF 1: {degrees_of_freedom(self)}") # condenser mode with cooling water estimation cold_side_outlet_temperature = self.cold_side_outlet.temperature[0].value self.hot_side_inlet.fix() self.cold_side_inlet.fix() - self.cold_side_outlet.unfix() + #self.hot_side_outlet.unfix() + self.cold_side_outlet.temperature.unfix() + # condenser mode without cooling water estimation + self.outlet_liquid_mass_balance.deactivate() + self.outlet_pressure_sat.deactivate() + print(f"DOF 2: {degrees_of_freedom(self)}") super().initialize_build(*args, **kwargs) self.area.unfix() self.cold_side_outlet.temperature[0].fix(cold_side_outlet_temperature) - self.cold_side.properties_in[0].mass_frac_phase_comp["Liq", "TDS"] - self.cold_side.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].fix() + """ + for j in self.cold_side.config.property_package.solute_set: + self.cold_side.properties_in[0].mass_frac_phase_comp["Liq", j] + self.cold_side.properties_in[0].mass_frac_phase_comp["Liq", j].fix() + """ - for j in self.cold_side.config.property_package.component_list: - self.cold_side.properties_in[0].flow_mass_phase_comp["Liq", j].unfix() + #for j in self.cold_side.config.property_package.component_list: + self.cold_side.properties_in[0].flow_mass_phase_comp["Liq", "H2O"].unfix() self.outlet_liquid_mass_balance.activate() self.outlet_pressure_sat.activate() @@ -177,19 +200,21 @@ def initialize_build(self, *args, **kwargs): ) ) - self.cold_side.properties_in[0].mass_frac_phase_comp["Liq", "TDS"].unfix() - self.cold_side.properties_in[0].flow_mass_phase_comp["Liq", "TDS"].fix() - - opt = get_solver(solver, optarg) - - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = opt.solve(self, tee=slc.tee) - init_log.info( - "Initialization Complete (w/ cooling water estimation): {}".format( - idaeslog.condition(res) + #for j in self.cold_side.config.property_package.solute_set: + # self.cold_side.properties_in[0].mass_frac_phase_comp["Liq", j].unfix() + # self.cold_side.properties_in[0].flow_mass_phase_comp["Liq", j].fix() + print(f"DOF 3: {degrees_of_freedom(self)}") + """ + opt = get_solver(solver, optarg) + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + init_log.info( + "Initialization Complete (w/ cooling water estimation): {}".format( + idaeslog.condition(res) + ) ) - ) - + """ @property def default_costing_method(self): return cost_heat_exchanger