diff --git a/watertap/flowsheets/RO_with_energy_recovery/RO_with_energy_recovery.py b/watertap/flowsheets/RO_with_energy_recovery/RO_with_energy_recovery.py index 3366297b8f..266ee2a770 100644 --- a/watertap/flowsheets/RO_with_energy_recovery/RO_with_energy_recovery.py +++ b/watertap/flowsheets/RO_with_energy_recovery/RO_with_energy_recovery.py @@ -23,8 +23,8 @@ from pyomo.network import Arc from idaes.core import FlowsheetBlock from watertap.core.solvers import get_solver +from idaes.core.util.initialization import propagate_state from idaes.core.util.model_statistics import degrees_of_freedom -from idaes.core.util.initialization import solve_indexed_blocks, propagate_state from idaes.models.unit_models import Mixer, Separator, Product, Feed from idaes.models.unit_models.mixer import MomentumMixingType from idaes.core import UnitModelCostingBlock @@ -42,6 +42,7 @@ from watertap.unit_models.pressure_exchanger import PressureExchanger from watertap.unit_models.pressure_changer import Pump, EnergyRecoveryDevice from watertap.core.util.initialization import assert_degrees_of_freedom +from watertap.tools.unit_models import calculate_operating_pressure from watertap.costing import WaterTAPCosting @@ -201,7 +202,7 @@ def build(erd_type=ERDtype.pressure_exchanger): def set_operating_conditions( m, water_recovery=0.5, - over_pressure=0.3, + over_pressure_factor=1.3, flow_vol=1e-3, salt_mass_conc=35e-3, solver=None, @@ -253,10 +254,10 @@ def set_operating_conditions( # pump 1, high pressure pump, 2 degrees of freedom (efficiency and outlet pressure) m.fs.P1.efficiency_pump.fix(0.80) # pump efficiency [-] operating_pressure = calculate_operating_pressure( - feed_state_block=m.fs.feed.properties[0], - over_pressure=over_pressure, - water_recovery=water_recovery, - NaCl_passage=0.01, + state_block=m.fs.feed.properties[0], + over_pressure_factor=over_pressure_factor, + water_recovery_mass=water_recovery, + salt_passage=0.01, solver=solver, ) m.fs.P1.control_volume.properties_out[0].pressure.fix(operating_pressure) @@ -312,53 +313,6 @@ def set_operating_conditions( ) -def calculate_operating_pressure( - feed_state_block=None, - over_pressure=0.15, - water_recovery=0.5, - NaCl_passage=0.01, - solver=None, -): - """ - Estimate operating pressure for RO unit model given the following arguments: - - Arguments: - feed_state_block: the state block of the RO feed that has the non-pressure state - variables initialized to their values (default=None) - over_pressure: the amount of operating pressure above the brine osmotic pressure - represented as a fraction (default=0.15) - water_recovery: the mass-based fraction of inlet H2O that becomes permeate - (default=0.5) - NaCl_passage: the mass-based fraction of inlet NaCl that becomes permeate - (default=0.01) - solver: solver object to be used (default=None) - """ - t = ConcreteModel() # create temporary model - prop = feed_state_block.config.parameters - t.brine = prop.build_state_block([0]) - - # specify state block - t.brine[0].flow_mass_phase_comp["Liq", "H2O"].fix( - value(feed_state_block.flow_mass_phase_comp["Liq", "H2O"]) - * (1 - water_recovery) - ) - t.brine[0].flow_mass_phase_comp["Liq", "NaCl"].fix( - value(feed_state_block.flow_mass_phase_comp["Liq", "NaCl"]) * (1 - NaCl_passage) - ) - t.brine[0].pressure.fix( - 101325 - ) # valid when osmotic pressure is independent of hydraulic pressure - t.brine[0].temperature.fix(value(feed_state_block.temperature)) - # calculate osmotic pressure - # since properties are created on demand, we must touch the property to create it - t.brine[0].pressure_osm_phase - # solve state block - results = solve_indexed_blocks(solver, [t.brine]) - assert_optimal_termination(results) - - return value(t.brine[0].pressure_osm_phase["Liq"]) * (1 + over_pressure) - - def solve(blk, solver=None, tee=False, check_termination=True): if solver is None: solver = get_solver() diff --git a/watertap/flowsheets/RO_with_energy_recovery/monte_carlo_sampling_RO_ERD.py b/watertap/flowsheets/RO_with_energy_recovery/monte_carlo_sampling_RO_ERD.py index 87a3051d0d..9f4b9a8f9c 100644 --- a/watertap/flowsheets/RO_with_energy_recovery/monte_carlo_sampling_RO_ERD.py +++ b/watertap/flowsheets/RO_with_energy_recovery/monte_carlo_sampling_RO_ERD.py @@ -69,7 +69,9 @@ def build_model( # Build, set, and initialize the system (these steps will change depending on the underlying model) m = build() - set_operating_conditions(m, water_recovery=0.5, over_pressure=0.3, solver=solver) + set_operating_conditions( + m, water_recovery=0.5, over_pressure_factor=1.3, solver=solver + ) initialize_system(m, solver=solver) # Check if we need to read in the default model values from a file diff --git a/watertap/flowsheets/RO_with_energy_recovery/tests/test_RO_with_energy_recovery_simulation.py b/watertap/flowsheets/RO_with_energy_recovery/tests/test_RO_with_energy_recovery_simulation.py index 936aa9078c..4320ff9e3b 100644 --- a/watertap/flowsheets/RO_with_energy_recovery/tests/test_RO_with_energy_recovery_simulation.py +++ b/watertap/flowsheets/RO_with_energy_recovery/tests/test_RO_with_energy_recovery_simulation.py @@ -142,7 +142,7 @@ def test_set_operating_conditions(self, system_frame): m = system_frame set_operating_conditions( - m, water_recovery=0.5, over_pressure=0.3, solver=solver + m, water_recovery=0.5, over_pressure_factor=1.3, solver=solver ) # check fixed variables @@ -513,7 +513,7 @@ def test_set_operating_conditions(self, system_frame): m = system_frame set_operating_conditions( - m, water_recovery=0.5, over_pressure=0.3, solver=solver + m, water_recovery=0.5, over_pressure_factor=1.3, solver=solver ) # check fixed variables diff --git a/watertap/tools/unit_models/__init__.py b/watertap/tools/unit_models/__init__.py new file mode 100644 index 0000000000..649a45a002 --- /dev/null +++ b/watertap/tools/unit_models/__init__.py @@ -0,0 +1,13 @@ +################################################################################# +# 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 .reverse_osmosis import * diff --git a/watertap/tools/unit_models/reverse_osmosis.py b/watertap/tools/unit_models/reverse_osmosis.py new file mode 100644 index 0000000000..298ec85b58 --- /dev/null +++ b/watertap/tools/unit_models/reverse_osmosis.py @@ -0,0 +1,116 @@ +################################################################################# +# 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/" +################################################################################# +""" +This module contains functions to be used with WaterTAP +ReverseOsmosis0D or ReverseOsmosis1D unit models. +""" + +from pyomo.environ import ( + ConcreteModel, + check_optimal_termination, + value, +) +from idaes.core.util.initialization import solve_indexed_blocks + +from watertap.property_models.seawater_prop_pack import SeawaterStateBlockData +from watertap.property_models.NaCl_prop_pack import NaClStateBlockData +from watertap.property_models.NaCl_T_dep_prop_pack import ( + NaClStateBlockData as NaClTDepStateBlockData, +) +from watertap.core import MembraneChannel0DBlock, MembraneChannel1DBlock +from watertap.core.solvers import get_solver + + +__all__ = ["calculate_operating_pressure"] + + +def calculate_operating_pressure( + state_block=None, + over_pressure_factor=1.15, + water_recovery_mass=0.5, + salt_passage=0, + solver=None, +): + """ + Estimate operating pressure for RO unit model given the following arguments: + + Arguments: + state_block: the state block of the RO feed that has the non-pressure state variables set to desired values (default=None) + over_pressure_factor: the amount of operating pressure above the brine osmotic pressure represented as a fraction (default=1.15) + water_recovery_mass: the mass-based fraction of inlet H2O that becomes permeate (default=0.5) + salt_passage: the mass-based fraction of inlet salt that becomes permeate (default=0) + solver: solver object to be used (default=None) + """ + + if any( + isinstance(state_block, cls) + for cls in [MembraneChannel0DBlock, MembraneChannel1DBlock] + ): + state_block = state_block.properties[0, 0] + + if not any( + isinstance(state_block, cls) + for cls in [SeawaterStateBlockData, NaClStateBlockData, NaClTDepStateBlockData] + ): + raise TypeError( + "state_block must be created with SeawaterParameterBlock, NaClParameterBlock, or NaClTDepParameterBlock" + ) + + if not 0 <= salt_passage < 0.999: + raise ValueError("salt_passage argument must be between 0 and 0.999") + + if not 1e-3 < water_recovery_mass < 0.999: + raise ValueError("water_recovery_mass argument must be between 0.001 and 0.999") + + if not over_pressure_factor >= 1.0: + raise ValueError( + "over_pressure_factor argument must be greater than or equal to 1.0" + ) + + comp = state_block.params.solute_set.first() + + if comp not in ["NaCl", "TDS"]: + raise ValueError( + f"salt_passage calculation only supported for NaCl or TDS components but found {comp}" + ) + + if solver is None: + solver = get_solver() + + tmp = ConcreteModel() # create temporary model + prop = state_block.config.parameters + + tmp.feed = prop.build_state_block([0]) + tmp.feed[0].pressure_osm_phase + + # specify state block + tmp.feed[0].flow_mass_phase_comp["Liq", "H2O"].fix( + value(state_block.flow_mass_phase_comp["Liq", "H2O"]) + * (1 - water_recovery_mass) + ) + tmp.feed[0].flow_mass_phase_comp["Liq", comp].fix( + value(state_block.flow_mass_phase_comp["Liq", comp]) * (1 - salt_passage) + ) + tmp.feed[0].temperature.fix(value(state_block.temperature)) + tmp.feed[0].pressure.fix(101325) + + # solve state block + results = solve_indexed_blocks(solver, [tmp.feed]) + + if not check_optimal_termination(results): + raise RuntimeError( + "Failed to solve temporary state block for operating pressure" + ) + + op_pressure = value(tmp.feed[0].pressure_osm_phase["Liq"]) * over_pressure_factor + + return op_pressure diff --git a/watertap/tools/unit_models/tests/test_reverse_osmosis_tools.py b/watertap/tools/unit_models/tests/test_reverse_osmosis_tools.py new file mode 100644 index 0000000000..40dc94d711 --- /dev/null +++ b/watertap/tools/unit_models/tests/test_reverse_osmosis_tools.py @@ -0,0 +1,220 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + +import pytest + +from pyomo.environ import ConcreteModel +from idaes.core import FlowsheetBlock +from idaes.models.unit_models import Feed + +from watertap.property_models.seawater_prop_pack import SeawaterParameterBlock +from watertap.property_models.NaCl_prop_pack import NaClParameterBlock +from watertap.property_models.NaCl_T_dep_prop_pack import ( + NaClParameterBlock as NaClTDepParameterBlock, +) +from watertap.property_models.multicomp_aq_sol_prop_pack import MCASParameterBlock + +from watertap.tools.unit_models import calculate_operating_pressure +from watertap.unit_models import ReverseOsmosis0D, ReverseOsmosis1D +from watertap.core.solvers import get_solver + +solver = get_solver() + + +def build_seawater_prop_model(): + """ + Create feed model using seawater property package + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + m.fs.feed = Feed(property_package=m.fs.properties) + m.fs.feed.properties[0].flow_mass_phase_comp["Liq", "H2O"].fix(0.965) + m.fs.feed.properties[0].flow_mass_phase_comp["Liq", "TDS"].fix(0.035) + m.fs.feed.properties[0].temperature.fix(273 + 25) + m.fs.feed.properties[0].pressure.fix(101325) + + m.fs.feed.initialize() + + return m + + +def build_nacl_prop_model(): + """ + Create feed model using NaCl property package + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = NaClParameterBlock() + m.fs.feed = Feed(property_package=m.fs.properties) + m.fs.feed.properties[0].flow_mass_phase_comp["Liq", "H2O"].fix(0.965) + m.fs.feed.properties[0].flow_mass_phase_comp["Liq", "NaCl"].fix(0.035) + m.fs.feed.properties[0].temperature.fix(273 + 25) + m.fs.feed.properties[0].pressure.fix(101325) + + m.fs.feed.initialize() + + return m + + +def build_nacl_tdep_prop_model(): + """ + Create feed model using NaCl temp dependence property package + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = NaClTDepParameterBlock() + m.fs.feed = Feed(property_package=m.fs.properties) + m.fs.feed.properties[0].flow_mass_phase_comp["Liq", "H2O"].fix(0.965) + m.fs.feed.properties[0].flow_mass_phase_comp["Liq", "NaCl"].fix(0.035) + m.fs.feed.properties[0].temperature.fix(273 + 25) + m.fs.feed.properties[0].pressure.fix(101325) + + m.fs.feed.initialize() + + return m + + +def build_ro0d_model(): + """ + Create RO0D model for testing + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = SeawaterParameterBlock() + + m.fs.RO = ReverseOsmosis0D( + property_package=m.fs.properties, + concentration_polarization_type="none", + mass_transfer_coefficient="none", + has_pressure_change=False, + ) + + m.fs.RO.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(0.965) + m.fs.RO.inlet.flow_mass_phase_comp[0, "Liq", "TDS"].fix(0.035) + m.fs.RO.inlet.temperature.fix(273 + 25) + m.fs.RO.inlet.pressure.fix(101325) + + return m + + +def build_ro1d_model(): + """ + Create RO1D model for testing + """ + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = NaClTDepParameterBlock() + + m.fs.RO = ReverseOsmosis1D( + property_package=m.fs.properties, + concentration_polarization_type="none", + mass_transfer_coefficient="none", + has_pressure_change=False, + ) + + m.fs.RO.inlet.flow_mass_phase_comp[0, "Liq", "H2O"].fix(0.965) + m.fs.RO.inlet.flow_mass_phase_comp[0, "Liq", "NaCl"].fix(0.035) + m.fs.RO.inlet.temperature.fix(273 + 25) + m.fs.RO.inlet.pressure.fix(101325) + + return m + + +@pytest.mark.component +def test_calculate_operating_pressure_sw(): + + m = build_seawater_prop_model() + osm = calculate_operating_pressure(m.fs.feed.properties[0]) + assert pytest.approx(osm, rel=1e-3) == 6098990.3 + + +@pytest.mark.component +def test_calculate_operating_pressure_nacl(): + + m = build_nacl_prop_model() + osm = calculate_operating_pressure(m.fs.feed.properties[0]) + assert pytest.approx(osm, rel=1e-3) == 6695261.0 + + +@pytest.mark.component +def test_calculate_operating_pressure_nacl_tdep(): + + m = build_nacl_tdep_prop_model() + osm = calculate_operating_pressure(m.fs.feed.properties[0]) + assert pytest.approx(osm, rel=1e-3) == 6677498.9 + + +@pytest.mark.component +def test_calculate_operating_pressure_ro0d(): + + m = build_ro0d_model() + osm1 = calculate_operating_pressure(m.fs.RO.feed_side.properties[0, 0]) + osm2 = calculate_operating_pressure(m.fs.RO.feed_side) + assert osm1 == osm2 + assert pytest.approx(osm1, rel=1e-3) == 6098990.3 + + +@pytest.mark.component +def test_calculate_operating_pressure_ro1d(): + + m = build_ro1d_model() + osm1 = calculate_operating_pressure(m.fs.RO.feed_side.properties[0, 0]) + osm2 = calculate_operating_pressure(m.fs.RO.feed_side) + assert osm1 == osm2 + assert pytest.approx(osm1, rel=1e-3) == 6677498.9 + + +@pytest.mark.unit +def test_calculate_operating_pressure_errors(): + + with pytest.raises( + TypeError, + match="state_block must be created with SeawaterParameterBlock, NaClParameterBlock, or NaClTDepParameterBlock", + ): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.properties = MCASParameterBlock(solute_list=["Na_+"]) + m.fs.feed = Feed(property_package=m.fs.properties) + calculate_operating_pressure(state_block=m.fs.feed.properties[0]) + + with pytest.raises( + TypeError, + match="state_block must be created with SeawaterParameterBlock, NaClParameterBlock, or NaClTDepParameterBlock", + ): + m = build_ro0d_model() + calculate_operating_pressure(state_block=m.fs.RO.inlet) + + with pytest.raises( + ValueError, match="salt_passage argument must be between 0 and 0.999" + ): + m = build_seawater_prop_model() + calculate_operating_pressure( + state_block=m.fs.feed.properties[0], salt_passage=1.1 + ) + + with pytest.raises( + ValueError, match="water_recovery_mass argument must be between 0.001 and 0.999" + ): + m = build_nacl_prop_model() + calculate_operating_pressure( + state_block=m.fs.feed.properties[0], water_recovery_mass=2.5 + ) + + with pytest.raises( + ValueError, + match="over_pressure_factor argument must be greater than or equal to 1.0", + ): + m = build_nacl_prop_model() + calculate_operating_pressure( + state_block=m.fs.feed.properties[0], over_pressure_factor=0.9 + )