IX model small fix for anions in demo (#1420)

* abs value of charge in eq_mass_removed

* add anion test for ix demo

* remove value call inside abs

Co-authored-by: bknueven <30801372+bknueven@users.noreply.github.com>

* remove value from imports

* add regen block to IX demo

* add requires_idaes_solver; test less things

---------

Co-authored-by: bknueven <30801372+bknueven@users.noreply.github.com>

watertap/examples/flowsheets/ion_exchange/ion_exchange_demo.py

@@ -19,7 +19,6 @@
 from pyomo.network import Arc
 
 from idaes.core import FlowsheetBlock, UnitModelCostingBlock
-from watertap.core.solvers import get_solver
 from idaes.core.util.model_statistics import degrees_of_freedom
 from idaes.core.util.scaling import calculate_scaling_factors
 from idaes.core.util.initialization import propagate_state
@@ -29,6 +28,7 @@
 from watertap.property_models.multicomp_aq_sol_prop_pack import MCASParameterBlock
 from watertap.unit_models.ion_exchange_0D import IonExchange0D
 from watertap.costing import WaterTAPCosting
+from watertap.core.solvers import get_solver
 
 import math
 
@@ -106,6 +106,7 @@ def ix_build(ions, target_ion=None, hazardous_waste=False, regenerant="NaCl"):
     # The must use the same property package as the ion exchange model
     m.fs.feed = Feed(property_package=m.fs.properties)
     m.fs.product = Product(property_package=m.fs.properties)
+    m.fs.regen = Product(property_package=m.fs.properties)
 
     # Configuration dictionary used to instantiate the ion exchange model:
     #   "property_package" indicates which property package to use for the ion exchange model
@@ -122,9 +123,11 @@ def ix_build(ions, target_ion=None, hazardous_waste=False, regenerant="NaCl"):
     # Add the ion exchange model to the flowsheet
     m.fs.ion_exchange = ix = IonExchange0D(**ix_config)
 
-    # Touch concentration properties so they are available for reporting.
+    # Touch properties so they are available for scaling, initialization, and reporting.
     ix.process_flow.properties_in[0].conc_mass_phase_comp[...]
     ix.process_flow.properties_out[0].conc_mass_phase_comp[...]
+    ix.regeneration_stream[0].conc_mass_phase_comp[...]
+    m.fs.feed.properties[0].flow_vol_phase[...]
     m.fs.feed.properties[0].conc_mass_phase_comp[...]
     m.fs.product.properties[0].conc_mass_phase_comp[...]
 
@@ -150,6 +153,7 @@ def ix_build(ions, target_ion=None, hazardous_waste=False, regenerant="NaCl"):
     # For example, in this next line the outlet Port on the Feed model is connected to the inlet Port on the ion exchange model
     m.fs.feed_to_ix = Arc(source=m.fs.feed.outlet, destination=ix.inlet)
     m.fs.ix_to_product = Arc(source=ix.outlet, destination=m.fs.product.inlet)
+    m.fs.ix_to_regen = Arc(source=ix.regen, destination=m.fs.regen.inlet)
 
     TransformationFactory("network.expand_arcs").apply_to(m)
 
@@ -209,11 +213,13 @@ def initialize_system(m):
     propagate_state(m.fs.feed_to_ix)
     # ... and then initialize the ion exchange model.
     m.fs.ion_exchange.initialize()
-    # With the ion exchange model initialized, we have initial guesses for the Product block
-    # and can propagate the state of the IX effluent stream.
+    # With the ion exchange model initialized, we have initial guesses for the Product and Regen blocks
+    # and can propagate the state of the IX effluent and regeneration stream.
     propagate_state(m.fs.ix_to_product)
-    # Finally, we initialize the product and costing blocks.
+    propagate_state(m.fs.ix_to_regen)
+    # Finally, we initialize the product, regen and costing blocks.
     m.fs.product.initialize()
+    m.fs.regen.initialize()
     m.fs.costing.initialize()
 
 
@@ -232,19 +238,22 @@ def optimize_system(m):
     # and (hopefully) a lower cost.
     ix.process_flow.properties_out[0].conc_mass_phase_comp["Liq", target_ion].fix(0.025)
 
-    # With the new effluent conditions for our ion exchange model, this will have implications for our downstream models (the Product block)
+    # With the new effluent conditions for our ion exchange model, this will have implications for our downstream models (the Product and Regen blocks)
     # Thus, we must re-propagate the new effluent state to these models...
     propagate_state(m.fs.ix_to_product)
+    propagate_state(m.fs.ix_to_regen)
     # ...and re-initialize them to our new conditions.
     m.fs.product.initialize()
+    m.fs.regen.initialize()
 
     # To adjust solution to fixed-pattern to achieve desired effluent, must unfix dimensionless_time.
     ix.dimensionless_time.unfix()
     # Can optimize around different design variables, e.g., bed_depth, service_flow_rate (or combinations of these)
     # Here demonstrates optimization around column design
     ix.number_columns.unfix()
     ix.bed_depth.unfix()
-    solver.solve(m)
+    optimized_results = solver.solve(m)
+    assert_optimal_termination(optimized_results)
 
 
 def get_ion_config(ions):

watertap/examples/flowsheets/ion_exchange/tests/test_ion_exchange_demo.py

@@ -40,23 +40,22 @@
 
 __author__ = "Kurban Sitterley"
 
-target_ion = "Ca_2+"
-ions = [target_ion]
-mass_frac = 1e-4
-feed_mass_frac = {target_ion: mass_frac}
+
 solver = get_solver()
 
 
-class TestIXDemo:
-    @pytest.fixture(scope="class")
-    def ix_0D(self):
+class TestIXDemoCa:
 
+    @pytest.fixture(scope="class")
+    def ix_0D_Ca(self):
+        target_ion = "Ca_2+"
+        ions = [target_ion]
         m = ixf.ix_build(ions)
         return m
 
     @pytest.mark.unit
-    def test_build_model(self, ix_0D):
-        m = ix_0D
+    def test_build_model(self, ix_0D_Ca):
+        m = ix_0D_Ca
 
         # Test basic build
         assert isinstance(m, ConcreteModel)
@@ -110,9 +109,9 @@ def test_build_model(self, ix_0D):
         assert_units_consistent(m)
 
     @pytest.mark.component
-    def test_specific_operating_conditions(self, ix_0D):
+    def test_specific_operating_conditions(self, ix_0D_Ca):
 
-        m = ix_0D
+        m = ix_0D_Ca
         ixf.set_operating_conditions(m)
         ixf.initialize_system(m)
         assert degrees_of_freedom(m) == 0
@@ -230,8 +229,8 @@ def test_specific_operating_conditions(self, ix_0D):
 
     @pytest.mark.component
     @pytest.mark.requires_idaes_solver
-    def test_optimization(self, ix_0D):
-        m = ix_0D
+    def test_optimization(self, ix_0D_Ca):
+        m = ix_0D_Ca
         ixf.optimize_system(m)
         isinstance(m.fs.obj, Objective)
         assert m.fs.obj.expr == m.fs.costing.LCOW
@@ -267,54 +266,25 @@ def test_optimization(self, ix_0D):
         assert degrees_of_freedom(m) == 0
 
         results_dict = {
-            "resin_diam": 0.0007,
-            "resin_bulk_dens": 0.7,
             "resin_surf_per_vol": 4285.714,
             "c_norm": {"Ca_2+": 0.99398},
-            "bed_vol_tot": 11.999,
             "bed_depth": 2.211,
-            "bed_porosity": 0.5,
-            "col_height": 4.2371,
-            "col_diam": 1.175,
-            "col_height_to_diam_ratio": 3.604409,
-            "number_columns": 5.0,
-            "t_breakthru": 133829.353,
-            "ebct": 239.999,
+            "t_breakthru": 133829.3,
             "vel_bed": 0.009213559,
-            "service_flow_rate": 15.0,
-            "N_Re": 6.449491,
+            "N_Re": 6.4494,
             "N_Sc": {"Ca_2+": 1086.956},
             "N_Sh": {"Ca_2+": 28.755},
             "N_Pe_particle": 0.122332,
             "N_Pe_bed": 386.44,
-            "resin_max_capacity": 3.0,
             "resin_eq_capacity": 2.9873,
             "resin_unused_capacity": 0.01266295,
-            "langmuir": {"Ca_2+": 0.7},
             "mass_removed": {"Ca_2+": 12546.81},
             "num_transfer_units": 38.872,
             "dimensionless_time": 1.3321,
             "partition_ratio": 418.227,
             "fluid_mass_transfer_coeff": {"Ca_2+": 3.7792e-05},
             "pressure_drop": 16.049,
             "bed_vol": 2.399,
-            "t_rinse": 1199.9,
-            "t_waste": 3600.0,
-            "t_cycle": 137429.353,
-            "regen_pump_power": 0.030195009,
-            "regen_tank_vol": 29.999,
-            "bw_flow": 0.0075372,
-            "bed_expansion_frac": 0.46395,
-            "rinse_flow": 0.049999999,
-            "bw_pump_power": 0.004551716,
-            "rinse_pump_power": 0.060390018,
-            "bed_expansion_h": 1.0259,
-            "main_pump_power": 6.7349,
-            "col_vol_per": 4.5988,
-            "col_vol_tot": 22.994,
-            "t_contact": 119.999,
-            "vel_inter": 0.018427119,
-            "bv_calc": 557.622,
             "lh": 12.909,
             "separation_factor": {"Ca_2+": 1.4285},
             "rate_coeff": {"Ca_2+": 0.0002313},
@@ -357,6 +327,7 @@ def test_optimization(self, ix_0D):
                 assert pytest.approx(r, rel=1e-3) == value(mv)
 
     @pytest.mark.unit
+    @pytest.mark.requires_idaes_solver
     def test_main_fun(self):
         m = ixf.main()
 
@@ -369,7 +340,6 @@ def test_main_fun(self):
             "bed_porosity": 0.5,
             "col_height": 4.2371,
             "col_diam": 1.1755,
-            "col_height_to_diam_ratio": 3.6044,
             "number_columns": 5.0,
             "t_breakthru": 133829.3,
             "ebct": 240,
@@ -380,7 +350,6 @@ def test_main_fun(self):
             "N_Sh": {"Ca_2+": 28.755},
             "N_Pe_particle": 0.122332,
             "N_Pe_bed": 386.4407,
-            "resin_max_capacity": 3.0,
             "resin_eq_capacity": 2.9873,
             "langmuir": {"Ca_2+": 0.7},
             "mass_removed": {"Ca_2+": 12546.81},
@@ -423,3 +392,163 @@ def test_main_fun(self):
                     assert pytest.approx(s, rel=1e-3) == value(mv[i])
             else:
                 assert pytest.approx(r, rel=1e-3) == value(mv)
+
+
+class TestIXDemoSO4:
+
+    @pytest.fixture(scope="class")
+    def ix_0D_SO4(self):
+        target_ion = "SO4_2-"
+        ions = [target_ion]
+        m = ixf.ix_build(ions)
+        return m
+
+    @pytest.mark.unit
+    def test_build_model(self, ix_0D_SO4):
+        m = ix_0D_SO4
+
+        # Test basic build
+        assert isinstance(m, ConcreteModel)
+        assert isinstance(m.fs, FlowsheetBlock)
+        assert isinstance(m.fs.properties, MCASParameterBlock)
+        assert isinstance(m.fs.costing, Block)
+        assert isinstance(m.fs.feed, Feed)
+        assert isinstance(m.fs.ion_exchange, IonExchange0D)
+        assert isinstance(m.fs.product, Product)
+
+        # Test port
+        assert isinstance(m.fs.feed.outlet, Port)
+        assert isinstance(m.fs.ion_exchange.inlet, Port)
+        assert isinstance(m.fs.ion_exchange.outlet, Port)
+        assert isinstance(m.fs.product.inlet, Port)
+
+        # # Test consting
+        assert isinstance(m.fs.ion_exchange.costing, Block)
+        assert isinstance(m.fs.ion_exchange.costing.capital_cost, Var)
+        assert isinstance(m.fs.ion_exchange.costing.fixed_operating_cost, Var)
+
+        var_str_list = [
+            "total_capital_cost",
+            "total_operating_cost",
+        ]
+        for var_str in var_str_list:
+            var = getattr(m.fs.costing, var_str)
+            assert isinstance(var, Var)
+
+        # Test arcs
+        arc_dict = {
+            m.fs.feed_to_ix: (m.fs.feed.outlet, m.fs.ion_exchange.inlet),
+            m.fs.ix_to_product: (m.fs.ion_exchange.outlet, m.fs.product.inlet),
+        }
+        for arc, port_tpl in arc_dict.items():
+            assert arc.source is port_tpl[0]
+            assert arc.destination is port_tpl[1]
+
+        # test configrations
+        assert len(m.fs.ion_exchange.config) == 11
+        assert not m.fs.ion_exchange.config.dynamic
+        assert not m.fs.ion_exchange.config.has_holdup
+        assert m.fs.ion_exchange.config.target_ion == "SO4_2-"
+        assert m.fs.ion_exchange.ion_exchange_type == IonExchangeType.anion
+        assert m.fs.ion_exchange.config.regenerant == RegenerantChem.NaCl
+        assert m.fs.ion_exchange.config.isotherm == IsothermType.langmuir
+
+        assert m.fs.ion_exchange.config.property_package is m.fs.properties
+        assert "H2O" in m.fs.properties.component_list
+
+        assert_units_consistent(m)
+
+    @pytest.mark.component
+    def test_specific_operating_conditions(self, ix_0D_SO4):
+
+        m = ix_0D_SO4
+        ixf.set_operating_conditions(m)
+        ixf.initialize_system(m)
+        assert degrees_of_freedom(m) == 0
+
+        solver = get_solver()
+        results = solver.solve(m)
+        assert_optimal_termination(results)
+        assert value(m.fs.feed.properties[0].flow_vol_phase["Liq"]) == pytest.approx(
+            0.05, rel=1e-3
+        )
+        assert value(m.fs.product.properties[0].flow_vol_phase["Liq"]) == pytest.approx(
+            0.049995010, rel=1e-3
+        )
+        assert value(
+            sum(
+                m.fs.feed.properties[0].conc_mass_phase_comp["Liq", j]
+                for j in m.fs.properties.ion_set
+            )
+        ) == pytest.approx(0.1, rel=1e-3)
+        assert value(
+            sum(
+                m.fs.product.properties[0].conc_mass_phase_comp["Liq", j]
+                for j in m.fs.properties.ion_set
+            )
+        ) == pytest.approx(8.821e-5, rel=1e-3)
+
+        results_dict = {
+            "resin_surf_per_vol": 4285.714,
+            "c_norm": {"SO4_2-": 0.473},
+            "bed_vol_tot": 12.0,
+            "col_height": 3.4887,
+            "col_diam": 1.498,
+            "col_height_to_diam_ratio": 2.327,
+            "number_columns": 4.0,
+            "t_breakthru": 136055.4,
+            "ebct": 240.0,
+            "vel_bed": 0.0070833,
+            "N_Re": 4.958333,
+            "N_Sc": {"SO4_2-": 943.3},
+            "N_Sh": {"SO4_2-": 25.095},
+            "N_Pe_particle": 0.107827,
+            "N_Pe_bed": 261.867,
+            "resin_max_capacity": 3.0,
+            "resin_eq_capacity": 1.685707,
+            "resin_unused_capacity": 1.314292,
+            "langmuir": {"SO4_2-": 0.7},
+            "mass_removed": {"SO4_2-": 7079.969},
+            "num_transfer_units": 39.087,
+            "dimensionless_time": 1.0,
+            "partition_ratio": 566.397,
+            "fluid_mass_transfer_coeff": {"SO4_2-": 3.8001e-05},
+            "pressure_drop": 9.450141,
+            "separation_factor": {"SO4_2-": 1.4285},
+            "rate_coeff": {"SO4_2-": 0.000232663},
+            "HTU": {"SO4_2-": 0.043492261},
+        }
+        for v, r in results_dict.items():
+            ixv = getattr(m.fs.ion_exchange, v)
+            if ixv.is_indexed():
+                for i, s in r.items():
+                    assert pytest.approx(s, rel=1e-3) == value(ixv[i])
+            else:
+                assert pytest.approx(r, rel=1e-3) == value(ixv)
+
+        sys_cost_results = {
+            "aggregate_capital_cost": 866570.3,
+            "aggregate_fixed_operating_cost": 5492.468,
+            "aggregate_variable_operating_cost": 0.0,
+            "aggregate_flow_electricity": 4.023,
+            "aggregate_flow_NaCl": 1016854.2,
+            "aggregate_flow_costs": {"electricity": 2468.7, "NaCl": 92576.0},
+            "total_capital_cost": 866570.3,
+            "total_operating_cost": 117029.8,
+            "aggregate_direct_capital_cost": 433285.1,
+            "maintenance_labor_chemical_operating_cost": 25997.1,
+            "total_fixed_operating_cost": 31489.578,
+            "total_variable_operating_cost": 85540.2,
+            "total_annualized_cost": 203686.8,
+            "annual_water_production": 1419950.1,
+            "LCOW": 0.14344,
+            "specific_energy_consumption": 0.022353,
+        }
+
+        for v, r in sys_cost_results.items():
+            mv = getattr(m.fs.costing, v)
+            if mv.is_indexed():
+                for i, s in r.items():
+                    assert pytest.approx(s, rel=1e-3) == value(mv[i])
+            else:
+                assert pytest.approx(r, rel=1e-3) == value(mv)

watertap/unit_models/ion_exchange_0D.py

@@ -1053,7 +1053,7 @@ def eq_partition_ratio(b):
                 self.target_ion_set, doc="Removed total mass of ion in equivalents"
             )
             def eq_mass_removed(b, j):
-                charge = prop_in.charge_comp[j]
+                charge = abs(prop_in.charge_comp[j])
                 return b.mass_removed[j] * charge == pyunits.convert(
                     b.resin_eq_capacity * b.resin_bulk_dens * b.bed_vol_tot,
                     to_units=pyunits.mol,