From f52db5e55771a98d20498b674f091d01a60463e3 Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Tue, 20 Aug 2024 14:48:09 -0400 Subject: [PATCH 1/9] add crystallizer flowsheet --- ...lizer_live_steam_with_condenser_chiller.py | 395 ++++++++++++++++++ 1 file changed, 395 insertions(+) create mode 100644 watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py 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..46ae065b01 --- /dev/null +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -0,0 +1,395 @@ +################################################################################# +# 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, +) +from pyomo.environ import ( + ConcreteModel, + Constraint, + Objective, + TransformationFactory, + 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 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_revised import Crystallization +from watertap.costing import WaterTAPCosting, CrystallizerCostType +from watertap.costing.unit_models.heat_exchanger import ( + cost_heat_exchanger, +) +from watertap.unit_models.steam_heater_0D import SteamHeater0D, Mode +from idaes.models.unit_models import Mixer, Product, Feed, Heater +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.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 idaes.models.unit_models.translator import Translator +from watertap.costing.unit_models.heater_chiller import ( + cost_heater_chiller, +) +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) + + #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) + + + # 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.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.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 + + 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.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=True + ) + + 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 + + + + # performance expressions + + # connections + m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.mixer.feed) + m.fs.s02 = Arc(source=m.fs.mixer.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.mixer.recycle) + + 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.condenser.hot_side_inlet) + m.fs.s09 = Arc(source=m.fs.condenser.hot_side_outlet, destination=m.fs.chiller.inlet) + + m.fs.eq_equal_temperature = Constraint( + expr=m.fs.chiller.control_volume.properties_out[0].temperature0 + == m.fs.condenser.cold_side.properties_in[0].temperature0 + ) + + + 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.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.mixer.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) + + m.fs.chiller.costing = UnitModelCostingBlock( + flowsheet_costing_block=m.fs.costing, + costing_method=cost_heater_chiller, + costing_method_arguments={"HC_type": "chiller"}, + ) + + m.fs.condenser.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) + + 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 + 50 ) + m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(0.1) + m.fs.heater.overall_heat_transfer_coefficient.fix(2e3) + m.fs.heater.area.fix(10) + + # 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(1) + 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 + 140) + 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.condenser.cold_side_inlet.pressure[0].fix(101325) + m.fs.condenser.cold_side_inlet.temperature[0].fix(273.15 + 10) + + + print("DOF finaleee:", 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() + + + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value + m.fs.mixer.recycle.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value + m.fs.mixer.recycle.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + + m.fs.mixer.initialize() + m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + propagate_state(m.fs.s02) + m.fs.heater.initialize() + m.fs.heater.cold_side_inlet.unfix() + 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() + propagate_state(m.fs.s07) + propagate_state(m.fs.s08) + m.fs.condenser.initialize() + propagate_state(m.fs.s09) + m.fs.chiller.initialize() + m.fs.costing.initialize() + + +def optimize_set_up(m): + # add objective + m.fs.objective = Objective(expr=m.fs.costing.LCOW) + + m.fs.heater.area.unfix() + m.fs.heater.area.setlb(1) + m.fs.heater.area.setub(750) + + # additional constraints + Temperature_rise = 10e5 + 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) + + #assert_degrees_of_freedom(m, 6) + + +def optimize(m, solver=None): + # --solve--- + return solve(m, solver=solver) + + +if __name__ == "__main__": + m = main() From 2b16d2e7f247fbe5296d923a1d040679b02e5cfc Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Wed, 21 Aug 2024 08:32:29 -0400 Subject: [PATCH 2/9] modifications --- ...lizer_live_steam_with_condenser_chiller.py | 58 +++++++++++++------ 1 file changed, 40 insertions(+), 18 deletions(-) diff --git a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py index 46ae065b01..ae630cfc4d 100644 --- a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -58,8 +58,8 @@ from watertap.costing.unit_models.heater_chiller import ( cost_heater_chiller, ) -from watertap.core.util.model_debug_mode import activate -activate() +#from watertap.core.util.model_debug_mode import activate +#activate() @@ -96,7 +96,7 @@ def build(): # Control volume flow blocks m.fs.feed = Feed(property_package= m.fs.properties_nacl) - m.fs.distillate = Product(property_package= m.fs.properties) + m.fs.distillate = Product(property_package= m.fs.properties_vapor) # unit models: steam heater, mixer, pump, crystalizer @@ -189,7 +189,7 @@ def eq_pressure(blk): mode=Mode.CONDENSER, estimate_cooling_water=True ) m.fs.chiller = Heater( - property_package=m.fs.properties_nacl, has_pressure_change=True + property_package=m.fs.properties_nacl, has_pressure_change=False ) m.fs.tb_vapor = Translator( @@ -197,12 +197,7 @@ def eq_pressure(blk): 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 ( @@ -233,11 +228,13 @@ def eq_pressure(blk): 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.condenser.hot_side_inlet) - m.fs.s09 = Arc(source=m.fs.condenser.hot_side_outlet, destination=m.fs.chiller.inlet) + + m.fs.s09 = Arc(source=m.fs.condenser.cold_side_outlet, destination=m.fs.chiller.inlet) + m.fs.s10 = Arc(source=m.fs.condenser.hot_side_outlet, destination=m.fs.distillate.inlet) m.fs.eq_equal_temperature = Constraint( - expr=m.fs.chiller.control_volume.properties_out[0].temperature0 - == m.fs.condenser.cold_side.properties_in[0].temperature0 + expr=m.fs.chiller.control_volume.properties_out[0].temperature + == m.fs.condenser.cold_side.properties_in[0].temperature ) @@ -264,6 +261,10 @@ def eq_pressure(blk): 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.calculate_scaling_factors(m) @@ -282,7 +283,7 @@ def add_costs(m): m.fs.chiller.costing = UnitModelCostingBlock( flowsheet_costing_block=m.fs.costing, - costing_method=cost_heater_chiller, + costing_method=cost_heater_chiller, costing_method_arguments={"HC_type": "chiller"}, ) @@ -321,11 +322,22 @@ def set_operating_conditions(m): 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.condenser.cold_side_inlet.pressure[0].fix(101325) m.fs.condenser.cold_side_inlet.temperature[0].fix(273.15 + 10) + m.fs.condenser.overall_heat_transfer_coefficient.fix(2e3) + m.fs.condenser.cold_side_inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].fix( + 0 + ) + m.fs.condenser.cold_side_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].set_value( + 12 + ) + m.fs.condenser.cold_side_outlet.temperature[0].fix(273.15 + 40) + m.fs.condenser.hot_side_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(0) + - print("DOF finaleee:", degrees_of_freedom(m.fs)) + print("DOF final:", degrees_of_freedom(m.fs)) def solve(blk, solver=None, tee=True): if solver is None: @@ -344,12 +356,11 @@ def initialize_system(m, solver=None, verbose=True): m.fs.feed.initialize() m.fs.crystallizer.initialize() - + propagate_state(m.fs.s01) m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value m.fs.mixer.recycle.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value m.fs.mixer.recycle.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value - m.fs.mixer.initialize() m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() propagate_state(m.fs.s02) @@ -364,12 +375,23 @@ def initialize_system(m, solver=None, verbose=True): m.fs.crystallizer.initialize() propagate_state(m.fs.s07) propagate_state(m.fs.s08) + 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.s09) m.fs.chiller.initialize() + propagate_state(m.fs.s10) + m.fs.distillate.initialize() m.fs.costing.initialize() + + def optimize_set_up(m): # add objective m.fs.objective = Objective(expr=m.fs.costing.LCOW) @@ -389,7 +411,7 @@ def optimize_set_up(m): def optimize(m, solver=None): # --solve--- return solve(m, solver=solver) - + if __name__ == "__main__": m = main() From bd85f8afdc9ac497ef91e4ac294eda93c4434d74 Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Thu, 22 Aug 2024 07:18:33 -0400 Subject: [PATCH 3/9] modifies costig --- watertap/costing/unit_models/crystallizer.py | 92 +---- .../crystallization/Crystallizer_MVR.py | 391 ++++++++++++++++++ 2 files changed, 392 insertions(+), 91 deletions(-) create mode 100644 watertap/flowsheets/crystallization/Crystallizer_MVR.py 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/flowsheets/crystallization/Crystallizer_MVR.py b/watertap/flowsheets/crystallization/Crystallizer_MVR.py new file mode 100644 index 0000000000..a8d4010591 --- /dev/null +++ b/watertap/flowsheets/crystallization/Crystallizer_MVR.py @@ -0,0 +1,391 @@ +################################################################################# +# 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 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_revised 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 +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) + + + # unit models: steam heater, mixer, pump, crystalizer, compressor + + 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.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.s02 = Arc(source=m.fs.mixer.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.mixer.recycle) + 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.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.inlet.temperature[0].fix(273.15 + 50 ) + m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(0.1) + m.fs.heater.overall_heat_transfer_coefficient.fix(2e3) + #m.fs.heater.area.fix(10) + + # 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(1) + + + 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.compressor.pressure_ratio.set_value(2) + m.fs.compressor.efficiency.fix(0.8) + + + + 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.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value + m.fs.mixer.recycle.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value + m.fs.mixer.recycle.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + m.fs.mixer.initialize() + m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + print(f"DOF: {degrees_of_freedom(m)}") + propagate_state(m.fs.s02) + propagate_state(m.fs.s09) + m.fs.heater.initialize() + #m.fs.heater.cold_side_inlet.unfix() + 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)) + m.fs.costing.initialize() + + + +def optimize_set_up(m): + # add objective + m.fs.objective = Objective(expr=m.fs.costing.LCOW) + + m.fs.crystallizer.inlet.temperature.unfix() + m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) + m.fs.crystallizer.inlet.temperature.setub(273.15 +55) + #m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].lb = 400 + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].set_value(800) + print("DOF final optimization:", degrees_of_freedom(m.fs)) + # additional constraints + Temperature_rise = 4 + m.fs.eq_heater_temperature_rise = Constraint( + expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] <= Temperature_rise) + + #assert_degrees_of_freedom(m, 6) + + +def optimize(m, solver=None): + # --solve--- + return solve(m, solver=solver) + + +if __name__ == "__main__": + m = main() + From deef8f9bc90977acc3a519df5a39b84acc539420 Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Thu, 22 Aug 2024 07:20:37 -0400 Subject: [PATCH 4/9] modifies steam heater --- watertap/unit_models/steam_heater_0D.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/watertap/unit_models/steam_heater_0D.py b/watertap/unit_models/steam_heater_0D.py index a6a5f9fe43..a2a206fb09 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__) @@ -154,12 +155,15 @@ def initialize_build(self, *args, **kwargs): self.hot_side_inlet.fix() self.cold_side_inlet.fix() self.cold_side_outlet.unfix() + self.area.unfix() + 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() @@ -177,9 +181,6 @@ 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: From f4d56c0e34ad8a3b109efabd9458da92df7af328 Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Thu, 22 Aug 2024 07:37:02 -0400 Subject: [PATCH 5/9] modifies chiller --- ...lizer_live_steam_with_condenser_chiller.py | 35 ++++++++++++------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py index ae630cfc4d..cadce2b9b0 100644 --- a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -72,8 +72,8 @@ def main(): set_operating_conditions(m) initialize_system(m, solver=solver) - optimize_set_up(m) - solve(m, solver=solver) + #optimize_set_up(m) + #solve(m, solver=solver) #print("\n***---optimization results---***") #display_system(m) @@ -230,13 +230,16 @@ def eq_pressure(blk): m.fs.s08 = Arc(source=m.fs.tb_vapor.outlet, destination=m.fs.condenser.hot_side_inlet) m.fs.s09 = Arc(source=m.fs.condenser.cold_side_outlet, destination=m.fs.chiller.inlet) + + m.fs.s10 = Arc(source=m.fs.condenser.hot_side_outlet, destination=m.fs.distillate.inlet) + m.fs.eq_equal_temperature = Constraint( - expr=m.fs.chiller.control_volume.properties_out[0].temperature - == m.fs.condenser.cold_side.properties_in[0].temperature + expr=m.fs.chiller.control_volume.heat[0] + == m.fs.condenser.cold_side.heat[0] ) - + TransformationFactory("network.expand_arcs").apply_to(m) add_costs(m) @@ -284,7 +287,7 @@ def add_costs(m): m.fs.chiller.costing = UnitModelCostingBlock( flowsheet_costing_block=m.fs.costing, costing_method=cost_heater_chiller, - costing_method_arguments={"HC_type": "chiller"}, + costing_method_arguments={"HC_type": "chiller"}, ) m.fs.condenser.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) @@ -303,10 +306,10 @@ def set_operating_conditions(m): m.fs.feed.temperature[0].fix(273.15 + 20) - m.fs.crystallizer.inlet.temperature[0].set_value(273.15 + 50 ) + m.fs.crystallizer.inlet.temperature[0].fix(273.15 + 100 ) m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(0.1) m.fs.heater.overall_heat_transfer_coefficient.fix(2e3) - m.fs.heater.area.fix(10) + m.fs.heater.area.set_value(10) # Fix m.fs.crystallizer.crystal_growth_rate.fix() @@ -366,6 +369,7 @@ def initialize_system(m, solver=None, verbose=True): propagate_state(m.fs.s02) m.fs.heater.initialize() m.fs.heater.cold_side_inlet.unfix() + m.fs.heater.area.unfix() 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 @@ -382,11 +386,14 @@ def initialize_system(m, solver=None, verbose=True): m.fs.condenser.initialize() m.fs.condenser.hot_side_inlet.unfix() + m.fs.condenser.hot_side_inlet.flow_mass_phase_comp[0,"Vap", "H2O"].fix() + propagate_state(m.fs.s09) - m.fs.chiller.initialize() + #m.fs.chiller.initialize() propagate_state(m.fs.s10) m.fs.distillate.initialize() + print("DOF final:", degrees_of_freedom(m.fs)) m.fs.costing.initialize() @@ -395,11 +402,13 @@ def initialize_system(m, solver=None, verbose=True): def optimize_set_up(m): # add objective m.fs.objective = Objective(expr=m.fs.costing.LCOW) - - m.fs.heater.area.unfix() - m.fs.heater.area.setlb(1) - m.fs.heater.area.setub(750) + m.fs.crystallizer.inlet.temperature.unfix() + m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) + m.fs.crystallizer.inlet.temperature.setub(273.15 +55) + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].lb = 400 + m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].set_value(800) + print("DOF final optimization:", degrees_of_freedom(m.fs)) # additional constraints Temperature_rise = 10e5 m.fs.eq_heater_temperature_rise = Constraint( From b790eabd556deeaa837f4be79662a03fb80ef96d Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Wed, 11 Sep 2024 11:05:37 -0400 Subject: [PATCH 6/9] add revised crystallizer model --- .../costing/unit_models/heat_exchanger.py | 3 +- .../crystallization/Crystallizer_MVR.py | 101 ++- ...lizer_live_steam_with_condenser_chiller.py | 245 +++--- watertap/unit_models/Crystallizer_revised.py | 831 ++++++++++++++++++ watertap/unit_models/steam_heater_0D.py | 54 +- 5 files changed, 1061 insertions(+), 173 deletions(-) create mode 100644 watertap/unit_models/Crystallizer_revised.py 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 index a8d4010591..2b3daead6e 100644 --- a/watertap/flowsheets/crystallization/Crystallizer_MVR.py +++ b/watertap/flowsheets/crystallization/Crystallizer_MVR.py @@ -11,17 +11,14 @@ ################################################################################# 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 @@ -31,7 +28,10 @@ from idaes.core.util.initialization import ( propagate_state, ) -from pyomo.util.check_units import assert_units_consistent +from watertap.unit_models.pressure_changer import Pump +from watertap.costing.unit_models.pump import ( + cost_pump, +) from idaes.core import FlowsheetBlock from idaes.core.util.model_statistics import degrees_of_freedom @@ -44,13 +44,8 @@ from watertap.property_models.unit_specific import cryst_prop_pack as props from watertap.unit_models.Crystallizer_revised 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 import 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 ( @@ -63,9 +58,7 @@ 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() @@ -80,15 +73,12 @@ def main(): set_operating_conditions(m) initialize_system(m, solver=solver) - #optimize_set_up(m) + optimize_set_up(m) - #interval_initializer(m) - #solve(m, solver=solver) + solve(m, solver=solver) - #print("\n***---optimization results---***") - #display_system(m) - #display_design(m) - #display_state(m) + display_system(m) + return m @@ -134,6 +124,8 @@ def build(): 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.tb_vapor = Translator( inlet_property_package=m.fs.properties, outlet_property_package=m.fs.properties_vapor, @@ -211,12 +203,16 @@ def eq_temperature(blk): @m.fs.tb_recycle.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.feed.flow_mass_phase_comp[0,"Liq", "NaCl"] == m.fs.crystallizer.solids.flow_mass_phase_comp[0,"Sol", "NaCl"] ) # connections m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.mixer.feed) - m.fs.s02 = Arc(source=m.fs.mixer.outlet, destination=m.fs.heater.cold_side_inlet) + 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) @@ -268,6 +264,11 @@ def add_costs(m): 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"}, ) @@ -285,26 +286,29 @@ def set_operating_conditions(m): m.fs.feed.temperature[0].fix(273.15 + 20) - m.fs.crystallizer.inlet.temperature[0].fix(273.15 + 50 ) - m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(0.1) + m.fs.crystallizer.temperature_operating.set_value(273.15 + 50 ) m.fs.heater.overall_heat_transfer_coefficient.fix(2e3) - #m.fs.heater.area.fix(10) + # 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(1) 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.temperature[0].fix(273.15 + 100) + m.fs.heater.cold_side_outlet.temperature[0].set_value(273.15 + 200) m.fs.crystallizer.inlet.pressure[0].set_value(101325) - m.fs.compressor.pressure_ratio.set_value(2) + m.fs.compressor.pressure_ratio.set_value(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) @@ -344,10 +348,17 @@ def initialize_system(m, solver=None, verbose=True): 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.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 @@ -359,6 +370,13 @@ def initialize_system(m, solver=None, verbose=True): 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() @@ -366,25 +384,34 @@ def initialize_system(m, solver=None, verbose=True): def optimize_set_up(m): # add objective m.fs.objective = Objective(expr=m.fs.costing.LCOW) + + m.fs.crystallizer.temperature_operating.unfix() + m.fs.crystallizer.temperature_operating.setlb(273.15 +30) + m.fs.crystallizer.temperature_operating.setub(273.15 +120) + + 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) + - m.fs.crystallizer.inlet.temperature.unfix() - m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) - m.fs.crystallizer.inlet.temperature.setub(273.15 +55) - #m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].lb = 400 - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].set_value(800) - print("DOF final optimization:", degrees_of_freedom(m.fs)) # additional constraints Temperature_rise = 4 m.fs.eq_heater_temperature_rise = Constraint( - expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] <= Temperature_rise) - - #assert_degrees_of_freedom(m, 6) + expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] <= 5) def optimize(m, solver=None): # --solve--- return solve(m, solver=solver) +def display_system(m): + print("Inlet temperature:", m.fs.crystallizer.inlet.temperature[0].value) + print("operating temperature:", m.fs.crystallizer.temperature_operating.value) + print("pressure ratio:", m.fs.compressor.pressure_ratio.value) + print("recycle H2O:", 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("Inlet temperature heater:", m.fs.heater.cold_side_inlet.temperature[0].value) + print("outlet temperature heater:", m.fs.heater.cold_side_outlet.temperature[0].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 index cadce2b9b0..40ca2d1b17 100644 --- a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -11,21 +11,27 @@ ################################################################################# 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 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 @@ -41,8 +47,13 @@ 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 Mixer, Product, Feed, Heater +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 ( @@ -54,10 +65,9 @@ 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 watertap.costing.unit_models.heater_chiller import ( - cost_heater_chiller, -) +import logging #from watertap.core.util.model_debug_mode import activate #activate() @@ -72,11 +82,13 @@ def main(): set_operating_conditions(m) initialize_system(m, solver=solver) - #optimize_set_up(m) - #solve(m, solver=solver) + optimize_set_up(m) + + #interval_initializer(m) + solve(m, solver=solver) #print("\n***---optimization results---***") - #display_system(m) + display_system(m) #display_design(m) #display_state(m) @@ -97,6 +109,10 @@ def build(): # Control volume flow blocks m.fs.feed = Feed(property_package= m.fs.properties_nacl) m.fs.distillate = Product(property_package= m.fs.properties_vapor) + m.fs.solids = Product(property_package=m.fs.properties) + + + # overall process water recovery (evaporation rate/solid content) # unit models: steam heater, mixer, pump, crystalizer @@ -128,6 +144,32 @@ def build(): 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["Liq", "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["Vap", "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 ( @@ -135,7 +177,7 @@ def eq_flow_mass_comp(blk): == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] ) @m.fs.tb_recycle.Constraint() - def eq_flow_mass_comp_vapor(blk): + 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"] @@ -161,7 +203,7 @@ def eq_flow_mass_comp(blk): == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] ) @m.fs.tb_heater_to_cryst.Constraint() - def eq_flow_mass_comp_vapor(blk): + 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"] @@ -175,45 +217,12 @@ def eq_temperature(blk): def eq_pressure(blk): return blk.properties_in[0].pressure == blk.properties_out[0].pressure - 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.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_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.eq_heater_temperature_rise = Constraint( + expr=m.fs.feed.flow_mass_phase_comp[0,"Liq", "NaCl"] == m.fs.crystallizer.solids.flow_mass_phase_comp[0,"Sol", "NaCl"] ) + + - @m.fs.tb_vapor.Constraint() - def eq_pressure(blk): - return blk.properties_in[0].pressure == blk.properties_out[0].pressure - - # performance expressions @@ -227,20 +236,13 @@ def eq_pressure(blk): m.fs.s06 = Arc(source=m.fs.tb_recycle.outlet, destination=m.fs.mixer.recycle) 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.condenser.hot_side_inlet) - - m.fs.s09 = Arc(source=m.fs.condenser.cold_side_outlet, destination=m.fs.chiller.inlet) - - - m.fs.s10 = Arc(source=m.fs.condenser.hot_side_outlet, destination=m.fs.distillate.inlet) + m.fs.s08 = Arc(source=m.fs.tb_vapor.outlet, destination=m.fs.distillate.inlet) + m.fs.s09 = Arc(source=m.fs.crystallizer.solids, destination=m.fs.solids.inlet) - - m.fs.eq_equal_temperature = Constraint( - expr=m.fs.chiller.control_volume.heat[0] - == m.fs.condenser.cold_side.heat[0] - ) - + #m.fs.eq_distillate = Constraint( + # expr=m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"] == m.fs.distillate.flow_mass_phase_comp[0,"Liq", "H2O"] ) + #m.fs.distillate.flow_mass_phase_comp[0,"Vap", "H2O"].fix(0) TransformationFactory("network.expand_arcs").apply_to(m) add_costs(m) @@ -264,10 +266,6 @@ def eq_pressure(blk): 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.calculate_scaling_factors(m) @@ -277,24 +275,18 @@ def add_costs(m): m.fs.costing = WaterTAPCosting() m.fs.crystallizer.costing = UnitModelCostingBlock( - flowsheet_costing_block=m.fs.costing, + 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, +) + 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.mixer.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) - m.fs.chiller.costing = UnitModelCostingBlock( - flowsheet_costing_block=m.fs.costing, - costing_method=cost_heater_chiller, - costing_method_arguments={"HC_type": "chiller"}, - ) - - m.fs.condenser.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) - 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): @@ -306,41 +298,39 @@ def set_operating_conditions(m): m.fs.feed.temperature[0].fix(273.15 + 20) - m.fs.crystallizer.inlet.temperature[0].fix(273.15 + 100 ) - m.fs.crystallizer.solids.flow_mass_phase_comp[0, "Sol", "NaCl"].fix(0.1) + 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(10) + 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(1) + 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 + 140) + m.fs.heater.hot_side_inlet.temperature.fix(273.15 + 70) + 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.temperature[0].set_value(273.15 + 80) m.fs.crystallizer.inlet.pressure[0].set_value(101325) - m.fs.condenser.cold_side_inlet.pressure[0].fix(101325) - m.fs.condenser.cold_side_inlet.temperature[0].fix(273.15 + 10) - m.fs.condenser.overall_heat_transfer_coefficient.fix(2e3) - m.fs.condenser.cold_side_inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].fix( - 0 - ) - m.fs.condenser.cold_side_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].set_value( - 12 - ) - m.fs.condenser.cold_side_outlet.temperature[0].fix(273.15 + 40) - m.fs.condenser.hot_side_inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(0) - + #m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) + #m.fs.crystallizer.inlet.temperature.setub(273.15 +55) + #m.fs.crystallizer.outlet.temperature.setlb(273.15 +50) + #m.fs.crystallizer.outlet.temperature.setub(273.15 +55) + m.fs.heater.area.setlb(10) + + - print("DOF final:", degrees_of_freedom(m.fs)) + print("DOF finaleee:", degrees_of_freedom(m.fs)) def solve(blk, solver=None, tee=True): if solver is None: @@ -359,60 +349,69 @@ def initialize_system(m, solver=None, verbose=True): m.fs.feed.initialize() m.fs.crystallizer.initialize() + + #propagate_state(m.fs.s07) + m.fs.distillate.initialize() + + propagate_state(m.fs.s01) + propagate_state(m.fs.s05) + propagate_state(m.fs.s06) m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value m.fs.mixer.recycle.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value m.fs.mixer.recycle.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + m.fs.mixer.initialize() m.fs.mixer.pressure_equality_constraints[0, 2].deactivate() + propagate_state(m.fs.s02) + m.fs.heater.initialize() m.fs.heater.cold_side_inlet.unfix() - m.fs.heater.area.unfix() + 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() + propagate_state(m.fs.s07) propagate_state(m.fs.s08) - 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() - m.fs.condenser.hot_side_inlet.flow_mass_phase_comp[0,"Vap", "H2O"].fix() - - - propagate_state(m.fs.s09) - #m.fs.chiller.initialize() - propagate_state(m.fs.s10) + m.fs.distillate.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"].value + m.fs.distillate.inlet.flow_mass_phase_comp[0,"Vap", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Liq", "H2O"].value + m.fs.distillate.inlet.temperature[0] = m.fs.crystallizer.vapor.temperature[0].value + m.fs.distillate.inlet.pressure[0] = m.fs.crystallizer.vapor.pressure[0].value m.fs.distillate.initialize() - print("DOF final:", degrees_of_freedom(m.fs)) - m.fs.costing.initialize() + propagate_state(m.fs.s09) + m.fs.solids.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 +55) + #m.fs.crystallizer.outlet.temperature.setlb(273.15 +50) + #m.fs.crystallizer.outlet.temperature.setub(273.15 +55) + 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) + + - m.fs.crystallizer.inlet.temperature.unfix() - m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) - m.fs.crystallizer.inlet.temperature.setub(273.15 +55) - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].lb = 400 - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"].set_value(800) - print("DOF final optimization:", degrees_of_freedom(m.fs)) # additional constraints - Temperature_rise = 10e5 + Temperature_rise = 4 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) + expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] <= 4) #assert_degrees_of_freedom(m, 6) @@ -420,7 +419,15 @@ def optimize_set_up(m): def optimize(m, solver=None): # --solve--- return solve(m, solver=solver) +def display_system(m): + print("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..49310b1352 --- /dev/null +++ b/watertap/unit_models/Crystallizer_revised.py @@ -0,0 +1,831 @@ +################################################################################# +# 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 + +_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.") + #self.inlet.temperature[0].unfix() + # --------------------------------------------------------------------- + # 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 + + state_args_out = deepcopy(state_args) + state_args_out["temperature"] = state_args["temperature"] + + self.properties_out.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args_out, + ) + + state_args_solids = deepcopy(state_args_out) + 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_out) + 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.") + # --------------------------------------------------------------------- + # 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/steam_heater_0D.py b/watertap/unit_models/steam_heater_0D.py index a2a206fb09..f9eb9991cd 100644 --- a/watertap/unit_models/steam_heater_0D.py +++ b/watertap/unit_models/steam_heater_0D.py @@ -94,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. @@ -109,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() @@ -140,33 +141,49 @@ 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() # 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.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.area.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() super().initialize_build(*args, **kwargs) self.area.unfix() self.cold_side_outlet.temperature[0].fix(cold_side_outlet_temperature) + """ 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() @@ -181,16 +198,21 @@ def initialize_build(self, *args, **kwargs): ) ) - 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 From 45d86d02ec40a6e04df17924be2ea2de975bd79f Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Thu, 12 Sep 2024 11:17:04 -0400 Subject: [PATCH 7/9] add purge to MVR config --- .../crystallization/Crystallizer_MVR.py | 100 +++++++++++------- ...lizer_live_steam_with_condenser_chiller.py | 35 ++---- watertap/unit_models/Crystallizer_revised.py | 20 ++-- watertap/unit_models/steam_heater_0D.py | 4 +- 4 files changed, 84 insertions(+), 75 deletions(-) diff --git a/watertap/flowsheets/crystallization/Crystallizer_MVR.py b/watertap/flowsheets/crystallization/Crystallizer_MVR.py index 2b3daead6e..c14a27498e 100644 --- a/watertap/flowsheets/crystallization/Crystallizer_MVR.py +++ b/watertap/flowsheets/crystallization/Crystallizer_MVR.py @@ -11,14 +11,17 @@ ################################################################################# 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 @@ -32,6 +35,7 @@ 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 @@ -44,8 +48,13 @@ from watertap.property_models.unit_specific import cryst_prop_pack as props from watertap.unit_models.Crystallizer_revised 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 Mixer, Product, Feed +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 ( @@ -58,9 +67,12 @@ 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 watertap.core.util.model_debug_mode import activate -activate() +from idaes.models.unit_models.separator import SplittingType +import logging +#from watertap.core.util.model_debug_mode import activate +#activate() @@ -78,7 +90,7 @@ def main(): solve(m, solver=solver) display_system(m) - + return m @@ -99,7 +111,7 @@ def build(): m.fs.distillate = Product(property_package= m.fs.properties_vapor) - # unit models: steam heater, mixer, pump, crystalizer, compressor + # unit models: steam heater, mixer, pump, crystalizer, compressor, separator m.fs.crystallizer = Crystallization(property_package=m.fs.properties) m.fs.heater = SteamHeater0D( @@ -126,6 +138,12 @@ def build(): 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, @@ -204,19 +222,18 @@ def eq_temperature(blk): 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.feed.flow_mass_phase_comp[0,"Liq", "NaCl"] == m.fs.crystallizer.solids.flow_mass_phase_comp[0,"Sol", "NaCl"] ) - - + # 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.mixer.recycle) + 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) @@ -287,8 +304,8 @@ def set_operating_conditions(m): 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() @@ -299,18 +316,21 @@ def set_operating_conditions(m): 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 + 100) - m.fs.heater.cold_side_outlet.temperature[0].set_value(273.15 + 200) + 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.set_value(1.5) + 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)) @@ -330,6 +350,7 @@ def initialize_system(m, solver=None, verbose=True): # 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 @@ -339,12 +360,21 @@ def initialize_system(m, solver=None, verbose=True): m.fs.compressor.initialize() propagate_state(m.fs.s01) + propagate_state(m.fs.s05) propagate_state(m.fs.s06) - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value - m.fs.mixer.recycle.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value - m.fs.mixer.recycle.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + + 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)}") @@ -382,37 +412,33 @@ def initialize_system(m, solver=None, verbose=True): def optimize_set_up(m): - # add objective 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.temperature_operating.unfix() - m.fs.crystallizer.temperature_operating.setlb(273.15 +30) - m.fs.crystallizer.temperature_operating.setub(273.15 +120) - + 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) - - # additional constraints - Temperature_rise = 4 - m.fs.eq_heater_temperature_rise = Constraint( - expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] <= 5) - def optimize(m, solver=None): # --solve--- return solve(m, solver=solver) def display_system(m): - print("Inlet temperature:", m.fs.crystallizer.inlet.temperature[0].value) print("operating temperature:", m.fs.crystallizer.temperature_operating.value) - print("pressure ratio:", m.fs.compressor.pressure_ratio.value) - print("recycle H2O:", 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("Inlet temperature heater:", m.fs.heater.cold_side_inlet.temperature[0].value) - print("outlet temperature heater:", m.fs.heater.cold_side_outlet.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 index 40ca2d1b17..1fcb152100 100644 --- a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -68,8 +68,8 @@ 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() +from watertap.core.util.model_debug_mode import activate +activate() @@ -219,9 +219,7 @@ def eq_pressure(blk): m.fs.eq_heater_temperature_rise = Constraint( - expr=m.fs.feed.flow_mass_phase_comp[0,"Liq", "NaCl"] == m.fs.crystallizer.solids.flow_mass_phase_comp[0,"Sol", "NaCl"] ) - - + expr=m.fs.heater.cold_side_outlet.temperature[0] - m.fs.heater.cold_side_inlet.temperature[0] == 4 * pyunits.K) # performance expressions @@ -311,7 +309,7 @@ def set_operating_conditions(m): 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 + 70) + 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) @@ -321,16 +319,9 @@ def set_operating_conditions(m): #m.fs.crystallizer.inlet.temperature[0].set_value(273.15 + 80) m.fs.crystallizer.inlet.pressure[0].set_value(101325) - #m.fs.crystallizer.inlet.temperature.setlb(273.15 +50) - #m.fs.crystallizer.inlet.temperature.setub(273.15 +55) - #m.fs.crystallizer.outlet.temperature.setlb(273.15 +50) - #m.fs.crystallizer.outlet.temperature.setub(273.15 +55) m.fs.heater.area.setlb(10) - - - - print("DOF finaleee:", degrees_of_freedom(m.fs)) + print("DOF final:", degrees_of_freedom(m.fs)) def solve(blk, solver=None, tee=True): if solver is None: @@ -401,31 +392,25 @@ def optimize_set_up(m): 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 +55) - #m.fs.crystallizer.outlet.temperature.setlb(273.15 +50) - #m.fs.crystallizer.outlet.temperature.setub(273.15 +55) + 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) - - - # additional constraints - Temperature_rise = 4 - 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) - #assert_degrees_of_freedom(m, 6) + 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("Inlet temperature:", m.fs.crystallizer.inlet.temperature[0].value) + 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("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) + diff --git a/watertap/unit_models/Crystallizer_revised.py b/watertap/unit_models/Crystallizer_revised.py index 49310b1352..580a3bbe33 100644 --- a/watertap/unit_models/Crystallizer_revised.py +++ b/watertap/unit_models/Crystallizer_revised.py @@ -41,6 +41,7 @@ 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__) @@ -486,8 +487,6 @@ 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 @@ -612,7 +611,6 @@ def initialize_build( hold_state=True, ) init_log.info_high("Initialization Step 1 Complete.") - #self.inlet.temperature[0].unfix() # --------------------------------------------------------------------- # Initialize other state blocks # Set state_args from inlet state @@ -630,22 +628,19 @@ def initialize_build( else: state_args[k] = state_dict[k].value - state_args_out = deepcopy(state_args) - state_args_out["temperature"] = state_args["temperature"] - self.properties_out.initialize( outlvl=outlvl, optarg=optarg, solver=solver, - state_args=state_args_out, + state_args=state_args, ) - state_args_solids = deepcopy(state_args_out) + 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] + ]["Liq", j] elif p == "Liq" or p == "Vap": state_args_solids["flow_mass_phase_comp"][p, j] = 1e-8 self.properties_solids.initialize( @@ -655,12 +650,12 @@ def initialize_build( state_args=state_args_solids, ) - state_args_vapor = deepcopy(state_args_out) + 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] + ]["Liq", j] elif p == "Liq" or p == "Sol": state_args_vapor["flow_mass_phase_comp"][p, j] = 1e-8 self.properties_vapor.initialize( @@ -670,6 +665,8 @@ def initialize_build( 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: @@ -682,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() diff --git a/watertap/unit_models/steam_heater_0D.py b/watertap/unit_models/steam_heater_0D.py index f9eb9991cd..9e56119e74 100644 --- a/watertap/unit_models/steam_heater_0D.py +++ b/watertap/unit_models/steam_heater_0D.py @@ -146,11 +146,13 @@ def initialize_build(self, *args, **kwargs): 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() @@ -171,7 +173,7 @@ def initialize_build(self, *args, **kwargs): # 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() From 23584545c517979a703641b147d0a3fa3517dea5 Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Tue, 1 Oct 2024 11:38:51 -0400 Subject: [PATCH 8/9] add chiller for crystallizer vapor condensation --- ...lizer_live_steam_with_condenser_chiller.py | 166 ++++++++++++++---- watertap/unit_models/steam_heater_0D.py | 4 +- 2 files changed, 130 insertions(+), 40 deletions(-) diff --git a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py index 1fcb152100..b6c3389fc0 100644 --- a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -54,6 +54,11 @@ ) 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 ( @@ -62,6 +67,9 @@ 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 @@ -109,7 +117,7 @@ def build(): # Control volume flow blocks m.fs.feed = Feed(property_package= m.fs.properties_nacl) m.fs.distillate = Product(property_package= m.fs.properties_vapor) - m.fs.solids = Product(property_package=m.fs.properties) + # overall process water recovery (evaporation rate/solid content) @@ -131,6 +139,25 @@ def build(): 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, @@ -139,6 +166,14 @@ def build(): ) 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, @@ -152,13 +187,13 @@ def build(): @m.fs.tb_vapor.Constraint() def eq_flow_mass_comp(blk): return ( - blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] + 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["Vap", "H2O"] + blk.properties_in[0].flow_mass_phase_comp["Liq", "H2O"] == blk.properties_out[0].flow_mass_phase_comp["Liq", "H2O"] ) @@ -220,27 +255,32 @@ def eq_pressure(blk): 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.mixer.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.mixer.recycle) + 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.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.distillate.inlet) - m.fs.s09 = Arc(source=m.fs.crystallizer.solids, destination=m.fs.solids.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) + - #m.fs.eq_distillate = Constraint( - # expr=m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"] == m.fs.distillate.flow_mass_phase_comp[0,"Liq", "H2O"] ) + - #m.fs.distillate.flow_mass_phase_comp[0,"Vap", "H2O"].fix(0) + TransformationFactory("network.expand_arcs").apply_to(m) add_costs(m) @@ -265,6 +305,12 @@ def eq_pressure(blk): 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 @@ -279,12 +325,26 @@ def add_costs(m): 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): @@ -319,7 +379,24 @@ def set_operating_conditions(m): #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.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)) @@ -341,28 +418,36 @@ def initialize_system(m, solver=None, verbose=True): m.fs.crystallizer.initialize() - #propagate_state(m.fs.s07) - m.fs.distillate.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.s05) - propagate_state(m.fs.s06) - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "H2O"].value - m.fs.mixer.recycle.flow_mass_phase_comp[0,"Liq", "NaCl"] = m.fs.crystallizer.outlet.flow_mass_phase_comp[0,"Liq", "NaCl"].value - m.fs.mixer.recycle.temperature[0] = m.fs.crystallizer.outlet.temperature[0].value - m.fs.mixer.recycle.pressure[0] = m.fs.crystallizer.outlet.pressure[0].value + 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.s02) + propagate_state(m.fs.s04) m.fs.heater.initialize() m.fs.heater.cold_side_inlet.unfix() - propagate_state(m.fs.s03) - propagate_state(m.fs.s04) + 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 @@ -370,15 +455,20 @@ def initialize_system(m, solver=None, verbose=True): m.fs.crystallizer.initialize() - propagate_state(m.fs.s07) - propagate_state(m.fs.s08) - m.fs.distillate.inlet.flow_mass_phase_comp[0,"Liq", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Vap", "H2O"].value - m.fs.distillate.inlet.flow_mass_phase_comp[0,"Vap", "H2O"] = m.fs.crystallizer.vapor.flow_mass_phase_comp[0,"Liq", "H2O"].value - m.fs.distillate.inlet.temperature[0] = m.fs.crystallizer.vapor.temperature[0].value - m.fs.distillate.inlet.pressure[0] = m.fs.crystallizer.vapor.pressure[0].value - m.fs.distillate.initialize() propagate_state(m.fs.s09) - m.fs.solids.initialize() + 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() @@ -391,7 +481,7 @@ def optimize_set_up(m): 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 +55) + 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) diff --git a/watertap/unit_models/steam_heater_0D.py b/watertap/unit_models/steam_heater_0D.py index 9e56119e74..8f27a0de78 100644 --- a/watertap/unit_models/steam_heater_0D.py +++ b/watertap/unit_models/steam_heater_0D.py @@ -171,8 +171,8 @@ def initialize_build(self, *args, **kwargs): #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() + self.outlet_liquid_mass_balance.deactivate() + self.outlet_pressure_sat.deactivate() print(f"DOF 2: {degrees_of_freedom(self)}") super().initialize_build(*args, **kwargs) From eb4acd699dc7182dcd6d7db785178fa70c774747 Mon Sep 17 00:00:00 2001 From: ElmiraShamlou Date: Tue, 1 Oct 2024 11:50:05 -0400 Subject: [PATCH 9/9] transfer revised crystallizer model to the original crystallizer model file --- .../crystallization/Crystallizer_MVR.py | 2 +- ...lizer_live_steam_with_condenser_chiller.py | 2 +- watertap/unit_models/crystallizer.py | 238 +++++++++--------- 3 files changed, 126 insertions(+), 116 deletions(-) diff --git a/watertap/flowsheets/crystallization/Crystallizer_MVR.py b/watertap/flowsheets/crystallization/Crystallizer_MVR.py index c14a27498e..f40a26a072 100644 --- a/watertap/flowsheets/crystallization/Crystallizer_MVR.py +++ b/watertap/flowsheets/crystallization/Crystallizer_MVR.py @@ -46,7 +46,7 @@ from idaes.core import UnitModelCostingBlock from watertap.property_models.unit_specific import cryst_prop_pack as props -from watertap.unit_models.Crystallizer_revised import Crystallization +from watertap.unit_models.crystallizer import Crystallization from watertap.costing import WaterTAPCosting, CrystallizerCostType from io import StringIO diff --git a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py index b6c3389fc0..b79e256695 100644 --- a/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py +++ b/watertap/flowsheets/crystallization/crystallizer_live_steam_with_condenser_chiller.py @@ -42,7 +42,7 @@ from idaes.core import UnitModelCostingBlock from watertap.property_models.unit_specific import cryst_prop_pack as props -from watertap.unit_models.Crystallizer_revised import Crystallization +from watertap.unit_models.crystallizer import Crystallization from watertap.costing import WaterTAPCosting, CrystallizerCostType from watertap.costing.unit_models.heat_exchanger import ( cost_heat_exchanger, 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)"