Merge pull request #1 from lxhowl/issue934dewatering

Issue934dewatering

watertap/examples/flowsheets/case_studies/BSM2/BSM2_flowsheet2.py

@@ -9,7 +9,7 @@
 # information, respectively. These files are also available online at the URL
 # "https://github.com/watertap-org/watertap/"
 #################################################################################
-__author__ = "Alejandro Garciadiego"
+__author__ = "Alejandro Garciadiego, Xinhong Liu"
 
 import pyomo.environ as pyo
 from pyomo.environ import (
@@ -200,17 +200,17 @@ def mass_transfer_R4(self, t):
 
     # Feed Water Conditions
     m.fs.FeedWater.flow_vol.fix(20648 * pyo.units.m**3 / pyo.units.day)
-    m.fs.FeedWater.temperature.fix(298.15 * pyo.units.K)
+    m.fs.FeedWater.temperature.fix(308.15 * pyo.units.K)
     m.fs.FeedWater.pressure.fix(1 * pyo.units.atm)
     m.fs.FeedWater.conc_mass_comp[0, "S_I"].fix(27 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "S_S"].fix(58 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "X_I"].fix(92 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "X_S"].fix(363 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "X_BH"].fix(50 * pyo.units.g / pyo.units.m**3)
-    m.fs.FeedWater.conc_mass_comp[0, "X_BA"].fix(1e-6 * pyo.units.g / pyo.units.m**3)
-    m.fs.FeedWater.conc_mass_comp[0, "X_P"].fix(1e-6 * pyo.units.g / pyo.units.m**3)
-    m.fs.FeedWater.conc_mass_comp[0, "S_O"].fix(1e-6 * pyo.units.g / pyo.units.m**3)
-    m.fs.FeedWater.conc_mass_comp[0, "S_NO"].fix(1e-6 * pyo.units.g / pyo.units.m**3)
+    m.fs.FeedWater.conc_mass_comp[0, "X_BA"].fix(0 * pyo.units.g / pyo.units.m**3)
+    m.fs.FeedWater.conc_mass_comp[0, "X_P"].fix(0 * pyo.units.g / pyo.units.m**3)
+    m.fs.FeedWater.conc_mass_comp[0, "S_O"].fix(0 * pyo.units.g / pyo.units.m**3)
+    m.fs.FeedWater.conc_mass_comp[0, "S_NO"].fix(0 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "S_NH"].fix(23 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "S_ND"].fix(5 * pyo.units.g / pyo.units.m**3)
     m.fs.FeedWater.conc_mass_comp[0, "X_ND"].fix(16 * pyo.units.g / pyo.units.m**3)
@@ -256,7 +256,7 @@ def mass_transfer_R4(self, t):
     m.fs.CL1.split_fraction[0, "effluent", "S_ALK"].fix(0.48956)
 
     # Sludge purge separator
-    m.fs.SP6.split_fraction[:, "recycle"].fix(0.975)
+    m.fs.SP6.split_fraction[:, "recycle"].fix(0.985)
 
     # Outlet pressure from recycle pump
     m.fs.P1.outlet.pressure.fix(101325)
@@ -325,7 +325,6 @@ def mass_transfer_R4(self, t):
     m.fs.RADM.liquid_outlet.temperature.fix(308.15)
 
     # Apply scaling
-
     m.fs.stream2adm = Arc(
         source=m.fs.RADM.liquid_outlet, destination=m.fs.adm_asm.inlet
     )
@@ -344,61 +343,79 @@ def mass_transfer_R4(self, t):
     pyo.TransformationFactory("network.expand_arcs").apply_to(m)
 
     iscale.calculate_scaling_factors(m.fs)
-    m.fs.set_default_scaling("conc_mass_comp", 1)
 
     # Initialize flowsheet
     # Apply sequential decomposition - 1 iteration should suffice
     seq = SequentialDecomposition()
-    seq.options.select_tear_method = "heuristic"
-    seq.options.tear_method = "Wegstein"
+    # seq.options.select_tear_method = "heuristic"
+    seq.options.tear_method = "Direct"
     seq.options.iterLim = 1
+    seq.options.tear_set = [m.fs.stream2, m.fs.stream10adm]
 
     G = seq.create_graph(m)
-
-    # # Uncomment this code to see tear set and initialization order
-    heuristic_tear_set = seq.tear_set_arcs(G, method="heuristic")
+    # Uncomment this code to see tear set and initialization order
     order = seq.calculation_order(G)
-    for o in heuristic_tear_set:
-        print(o.name)
+    print("Initialization Order")
     for o in order:
         print(o[0].name)
 
     # Initial guesses for flow into first reactor
-    tear_guesses = {
-        "flow_vol": {0: 0.26},
+    tear_guesses1 = {
+        "flow_vol": {0: 103531 / 24 / 3600},
         "conc_mass_comp": {
             (0, "S_I"): 0.028,
-            (0, "S_S"): 0.059,
-            (0, "X_I"): 0.094,
-            (0, "X_S"): 0.356,
-            (0, "X_BH"): 0.05,
-            (0, "X_BA"): 0.0009,
-            (0, "X_P"): 0.0006,
-            (0, "S_O"): 0.000017,
-            (0, "S_NO"): 0.000117,
-            (0, "S_NH"): 0.034,
-            (0, "S_ND"): 0.005,
-            (0, "X_ND"): 0.015,
+            (0, "S_S"): 0.012,
+            (0, "X_I"): 1.532,
+            (0, "X_S"): 0.069,
+            (0, "X_BH"): 2.233,
+            (0, "X_BA"): 0.167,
+            (0, "X_P"): 0.964,
+            (0, "S_O"): 0.0011,
+            (0, "S_NO"): 0.0073,
+            (0, "S_NH"): 0.0072,
+            (0, "S_ND"): 0.0016,
+            (0, "X_ND"): 0.0040,
         },
-        "alkalinity": {0: 7.69},
-        "temperature": {0: 298.15},
+        "alkalinity": {0: 0.0052},
+        "temperature": {0: 308.15},
+        "pressure": {0: 101325},
+    }
+
+    tear_guesses2 = {
+        "flow_vol": {0: 178 / 24 / 3600},
+        "conc_mass_comp": {
+            (0, "S_I"): 0.028,
+            (0, "S_S"): 0.048,
+            (0, "X_I"): 10.362,
+            (0, "X_S"): 20.375,
+            (0, "X_BH"): 10.210,
+            (0, "X_BA"): 0.553,
+            (0, "X_P"): 3.204,
+            (0, "S_O"): 0.00025,
+            (0, "S_NO"): 0.00169,
+            (0, "S_NH"): 0.0289,
+            (0, "S_ND"): 0.00468,
+            (0, "X_ND"): 0.906,
+        },
+        "alkalinity": {0: 0.00715},
+        "temperature": {0: 308.15},
         "pressure": {0: 101325},
     }
 
     # Pass the tear_guess to the SD tool
-    seq.set_guesses_for(m.fs.R1.inlet, tear_guesses)
+    seq.set_guesses_for(m.fs.R1.inlet, tear_guesses1)
+    seq.set_guesses_for(m.fs.asm_adm.inlet, tear_guesses2)
 
     def function(unit):
         unit.initialize(outlvl=idaeslog.INFO_HIGH)
 
     seq.run(m, function)
 
-    # solver = get_solver(options={"bound_push": 1e-8,"max_iter":0})
     solver = get_solver()
     results = solver.solve(m, tee=True)
-    print(degrees_of_freedom(m))
-    # pyo.assert_optimal_termination(results)
-    # m.display()
+
+    pyo.assert_optimal_termination(results)
+    m.display()
 
     print(large_residuals_set(m))
 

watertap/examples/flowsheets/case_studies/BSM2/bsm2_test2.py

@@ -21,7 +21,7 @@
 """
 
 # Some more information about this module
-__author__ = "Alejandro Garciadiego"
+__author__ = "Alejandro Garciadiego, Xinhong Liu"
 
 import pytest
 import pytest

watertap/examples/flowsheets/case_studies/anaerobic_digestor/ADM1_flowsheet.py

@@ -51,7 +51,7 @@ def build_flowsheet():
     m.fs.R1.inlet.flow_vol.fix(170 * pyo.units.m**3 / pyo.units.day)
     m.fs.R1.inlet.temperature.fix(308.15 * pyo.units.K)
     m.fs.R1.inlet.pressure.fix(1 * pyo.units.atm)
-    m.fs.R1.inlet.conc_mass_comp[0, "S_su"].fix(1 * pyo.units.mg / pyo.units.liter)
+    m.fs.R1.inlet.conc_mass_comp[0, "S_su"].fix(10 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "S_aa"].fix(1 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "S_fa"].fix(1 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "S_va"].fix(1 * pyo.units.mg / pyo.units.liter)
@@ -72,7 +72,7 @@ def build_flowsheet():
     m.fs.R1.inlet.conc_mass_comp[0, "X_ch"].fix(5000 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "X_pr"].fix(20000 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "X_li"].fix(5000 * pyo.units.mg / pyo.units.liter)
-    m.fs.R1.inlet.conc_mass_comp[0, "X_su"].fix(1 * pyo.units.mg / pyo.units.liter)
+    m.fs.R1.inlet.conc_mass_comp[0, "X_su"].fix(0.0 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "X_aa"].fix(10 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "X_fa"].fix(10 * pyo.units.mg / pyo.units.liter)
     m.fs.R1.inlet.conc_mass_comp[0, "X_c4"].fix(10 * pyo.units.mg / pyo.units.liter)
@@ -90,9 +90,9 @@ def build_flowsheet():
     m.fs.R1.liquid_outlet.temperature.fix(308.15 * pyo.units.K)
 
     # TO DO: Fix initialization
-    m.fs.R1.initialize(outlvl=idaeslog.INFO_HIGH, optarg={"bound_push": 1e-8})
+    m.fs.R1.initialize(outlvl=idaeslog.INFO_HIGH)
 
-    solver = get_solver(options={"bound_push": 1e-8})
+    solver = get_solver()
 
     results = solver.solve(m, tee=True)
 

watertap/examples/flowsheets/case_studies/anaerobic_digestor/tests/test_adm1_flowsheet.py

@@ -64,7 +64,7 @@ def test_results(self, model):
         )
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_su"]
-        ) == pytest.approx(1.193e-2, rel=1e-2)
+        ) == pytest.approx(1.195e-2, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_aa"]
         ) == pytest.approx(5.31e-3, rel=1e-2)
@@ -82,25 +82,25 @@ def test_results(self, model):
         ) == pytest.approx(0.01578, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_ac"]
-        ) == pytest.approx(0.17319, rel=1e-2)
+        ) == pytest.approx(0.19763, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_h2"]
-        ) == pytest.approx(2.35e-7, rel=1e-2)
+        ) == pytest.approx(2.36e-7, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_ch4"]
-        ) == pytest.approx(0.0541, rel=1e-2)
+        ) == pytest.approx(0.05508, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_IC"]
-        ) == pytest.approx(0.15436 * 12, rel=1e-2)
+        ) == pytest.approx(0.15268 * 12, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_IN"]
-        ) == pytest.approx(0.13022 * 14, rel=1e-2)
+        ) == pytest.approx(0.13023 * 14, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "S_I"]
         ) == pytest.approx(0.32869, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "X_c"]
-        ) == pytest.approx(0.3086, rel=1e-2)
+        ) == pytest.approx(0.30869, rel=1e-2)
         assert value(
             model.fs.R1.liquid_outlet.conc_mass_comp[0, "X_ch"]
         ) == pytest.approx(0.02794, rel=1e-2)

watertap/property_models/activated_sludge/asm1_properties.py

@@ -36,8 +36,7 @@
 import idaes.core.util.scaling as iscale
 
 # Some more information about this module
-__author__ = "Andrew Lee, Adam Atia"
-
+__author__ = "Andrew Lee, Adam Atia, Xinhong Liu"
 
 # Set up logger
 _log = idaeslog.getLogger(__name__)
@@ -283,7 +282,7 @@ def build(self):
         self.temperature = pyo.Var(
             domain=pyo.NonNegativeReals,
             initialize=298.15,
-            bounds=(293.15, 323.15),
+            bounds=(273.15, 323.15),
             doc="Temperature",
             units=pyo.units.K,
         )
@@ -361,7 +360,7 @@ def energy_density_expression(self):
 
         iscale.set_scaling_factor(self.flow_vol, 1e1)
         iscale.set_scaling_factor(self.temperature, 1e-1)
-        iscale.set_scaling_factor(self.pressure, 1e-3)
+        iscale.set_scaling_factor(self.pressure, 1e-6)
         iscale.set_scaling_factor(self.conc_mass_comp, 1e1)
 
     def get_material_flow_terms(self, p, j):

watertap/property_models/activated_sludge/asm1_reactions.py

@@ -38,7 +38,7 @@
 
 
 # Some more information about this module
-__author__ = "Andrew Lee"
+__author__ = "Andrew Lee, Xinhong Liu"
 
 
 # Set up logger
@@ -465,13 +465,10 @@ def rate_expression_rule(b, r):
                     # R7: Hydrolysis of entrapped organics
                     return b.reaction_rate[r] == pyo.units.convert(
                         b.params.k_h
-                        * (b.conc_mass_comp_ref["X_S"] / b.conc_mass_comp_ref["X_BH"])
+                        * b.conc_mass_comp_ref["X_S"]
                         / (
-                            b.params.K_X
-                            + (
-                                b.conc_mass_comp_ref["X_S"]
-                                / b.conc_mass_comp_ref["X_BH"]
-                            )
+                            b.params.K_X * b.conc_mass_comp_ref["X_BH"]
+                            + b.conc_mass_comp_ref["X_S"]
                         )
                         * (
                             (
@@ -493,7 +490,10 @@ def rate_expression_rule(b, r):
                     # R8: Hydrolysis of entrapped organic nitrogen
                     return b.reaction_rate[r] == (
                         b.reaction_rate["R7"]
-                        * (b.conc_mass_comp_ref["X_ND"] / b.conc_mass_comp_ref["X_S"])
+                        * (
+                            b.conc_mass_comp_ref["X_ND"] / b.conc_mass_comp_ref["X_S"]
+                            + 1e-10
+                        )
                     )
                 else:
                     raise BurntToast()

watertap/property_models/anaerobic_digestion/adm1_properties.py

@@ -36,7 +36,7 @@
 import idaes.core.util.scaling as iscale
 
 # Some more information about this module
-__author__ = "Alejandro Garciadiego, Adam Atia"
+__author__ = "Alejandro Garciadiego, Adam Atia, Xinhong Liu"
 # Using Andrew Lee's formulation of ASM1 as a template
 
 # Set up logger
@@ -279,14 +279,40 @@ def build(self):
         self.temperature = pyo.Var(
             domain=pyo.NonNegativeReals,
             initialize=298.15,
-            bounds=(298.15, 323.15),
+            bounds=(273.15, 323.15),
             doc="Temperature",
             units=pyo.units.K,
         )
+        Comp_dict = {
+            "S_su": 0.012,
+            "S_aa": 0.0053,
+            "S_fa": 0.099,
+            "S_va": 0.012,
+            "S_bu": 0.013,
+            "S_pro": 0.016,
+            "S_ac": 0.20,
+            "S_h2": 2.3e-7,
+            "S_ch4": 0.055,
+            "S_IC": 0.15 * 12,
+            "S_IN": 0.13 * 14,
+            "S_I": 0.33,
+            "X_c": 0.31,
+            "X_ch": 0.028,
+            "X_pr": 0.10,
+            "X_li": 0.029,
+            "X_su": 0.42,
+            "X_aa": 1.18,
+            "X_fa": 0.24,
+            "X_c4": 0.43,
+            "X_pro": 0.14,
+            "X_ac": 0.76,
+            "X_h2": 0.32,
+            "X_I": 25.6,
+        }
         self.conc_mass_comp = pyo.Var(
             self.params.solute_set,
             domain=pyo.NonNegativeReals,
-            initialize=0.001,
+            initialize=Comp_dict,
             doc="Component mass concentrations",
             units=pyo.units.kg / pyo.units.m**3,
         )
@@ -366,12 +392,13 @@ def energy_density_expression(self):
             rule=energy_density_expression, doc="Energy density term"
         )
 
-        iscale.set_scaling_factor(self.flow_vol, 1e1)
+        iscale.set_scaling_factor(self.flow_vol, 1e5)
         iscale.set_scaling_factor(self.temperature, 1e-1)
-        iscale.set_scaling_factor(self.pressure, 1e-3)
-        iscale.set_scaling_factor(self.conc_mass_comp, 1e1)
-        iscale.set_scaling_factor(self.anions, 1e1)
-        iscale.set_scaling_factor(self.cations, 1e1)
+        iscale.set_scaling_factor(self.pressure, 1e-6)
+        iscale.set_scaling_factor(self.conc_mass_comp, 1e2)
+        iscale.set_scaling_factor(self.conc_mass_comp["S_h2"], 1e5)
+        iscale.set_scaling_factor(self.anions, 1e2)
+        iscale.set_scaling_factor(self.cations, 1e2)
 
     def get_material_flow_terms(self, p, j):
         return self.material_flow_expression[j]