Create unit model tools directory and move calculate_operating_pressure

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()

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

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

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 *

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