From 7d3589008814463cffe8bd6ee79d8cc8f7c7c8b6 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Thu, 16 Mar 2023 16:00:58 -0400 Subject: [PATCH 01/21] add --- .../flowsheets/storage_tank/storage_tank.py | 340 ++++++++++++++++++ ...boxic_activated_sludge_process_global.yaml | 46 +++ 2 files changed, 386 insertions(+) create mode 100644 watertap/examples/flowsheets/storage_tank/storage_tank.py create mode 100644 watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml diff --git a/watertap/examples/flowsheets/storage_tank/storage_tank.py b/watertap/examples/flowsheets/storage_tank/storage_tank.py new file mode 100644 index 0000000000..9ae3a55eb4 --- /dev/null +++ b/watertap/examples/flowsheets/storage_tank/storage_tank.py @@ -0,0 +1,340 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2023, 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 os +import idaes.logger as idaeslog +from pyomo.environ import ( + ConcreteModel, + Expression, + value, + TransformationFactory, + units as pyunits, +) +from pyomo.network import Arc, SequentialDecomposition +from pyomo.util.check_units import assert_units_consistent + +from idaes.core import FlowsheetBlock +from idaes.core.solvers import get_solver +import idaes.core.util.scaling as iscale +from idaes.models.unit_models import Product +from idaes.core import UnitModelCostingBlock + +from watertap.core.wt_database import Database +import watertap.core.zero_order_properties as prop_ZO +from watertap.core.util.initialization import assert_degrees_of_freedom, check_solve +from watertap.unit_models.zero_order import ( + FeedZO, + StorageTankZO, +) +from watertap.core.zero_order_costing import ZeroOrderCosting + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +def main(): + m = build() + + set_operating_conditions(m) + assert_degrees_of_freedom(m, 0) + assert_units_consistent(m) + + initialize_system(m) + + results = solve(m, checkpoint="solve flowsheet after initializing system") + display_results(m) + + add_costing(m) + initialize_costing(m) + assert_degrees_of_freedom(m, 0) + assert_units_consistent(m) + + results = solve(m, checkpoint="solve flowsheet after costing") + + display_metrics_results(m) + display_additional_results(m) + return m, results + + +def build(): + m = ConcreteModel() + m.db = Database() + + m.fs = FlowsheetBlock(dynamic=False) + m.fs.prop = prop_ZO.WaterParameterBlock( + solute_list=["bod", "tss", "tkn", "phosphorus"] + ) + + # components + m.fs.feed = FeedZO(property_package=m.fs.prop) + m.fs.product_H2O = Product(property_package=m.fs.prop) + + m.fs.suboxicASM = StorageTankZO(property_package=m.fs.prop, database=m.db) + + # connections + m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.suboxicASM.inlet) + m.fs.s02 = Arc(source=m.fs.suboxicASM.outlet, destination=m.fs.product_H2O.inlet) + + TransformationFactory("network.expand_arcs").apply_to(m) + + # scale + iscale.calculate_scaling_factors(m) + return m + + +def set_operating_conditions(m): + # feed + flow_vol = 10 * pyunits.Mgallons / pyunits.day + conc_bod = 314 * pyunits.mg / pyunits.L + conc_tss = 336 * pyunits.mg / pyunits.L + conc_tkn = 52 * pyunits.mg / pyunits.L + conc_phosphorus = 6 * pyunits.mg / pyunits.L + + m.fs.feed.flow_vol[0].fix(flow_vol) + m.fs.feed.conc_mass_comp[0, "bod"].fix(conc_bod) + m.fs.feed.conc_mass_comp[0, "tss"].fix(conc_tss) + m.fs.feed.conc_mass_comp[0, "tkn"].fix(conc_tkn) + m.fs.feed.conc_mass_comp[0, "phosphorus"].fix(conc_phosphorus) + solve(m.fs.feed, checkpoint="solve feed block") + + # suboxicASM + m.fs.suboxicASM.load_parameters_from_database(use_default_removal=True) + + +def initialize_system(m): + seq = SequentialDecomposition() + seq.options.tear_set = [] + seq.options.iterLim = 1 + seq.run(m, lambda u: u.initialize()) + + +def solve(blk, solver=None, checkpoint=None, tee=False, fail_flag=True): + if solver is None: + solver = get_solver() + results = solver.solve(blk, tee=tee) + check_solve(results, checkpoint=checkpoint, logger=_log, fail_flag=fail_flag) + return results + + +def display_results(m): + unit_list = ["feed", "suboxicASM"] + for u in unit_list: + m.fs.component(u).report() + + +def add_costing(m): + source_file = os.path.join( + os.path.dirname(os.path.abspath(__file__)), + "suboxic_activated_sludge_process_global.yaml", + ) + m.fs.costing = ZeroOrderCosting(case_study_definition=source_file) + # typing aid + costing_kwargs = {"flowsheet_costing_block": m.fs.costing} + m.fs.suboxicASM.costing = UnitModelCostingBlock(**costing_kwargs) + + m.fs.costing.cost_process() + m.fs.costing.add_electricity_intensity(m.fs.product_H2O.properties[0].flow_vol) + m.fs.costing.add_LCOW(m.fs.product_H2O.properties[0].flow_vol) + + # other levelized costs + m.fs.costing.annual_water_inlet = Expression( + expr=m.fs.costing.utilization_factor + * pyunits.convert( + m.fs.feed.properties[0].flow_vol, + to_units=pyunits.m**3 / m.fs.costing.base_period, + ) + ) + + m.fs.costing.annual_water_production = Expression( + expr=m.fs.costing.utilization_factor + * pyunits.convert( + m.fs.product_H2O.properties[0].flow_vol, + to_units=pyunits.m**3 / m.fs.costing.base_period, + ) + ) + + m.fs.costing.total_annualized_cost = Expression( + expr=( + m.fs.costing.total_capital_cost * m.fs.costing.capital_recovery_factor + + m.fs.costing.total_operating_cost + ) + ) + + m.fs.costing.LCOT = Expression( + expr=(m.fs.costing.total_annualized_cost / m.fs.costing.annual_water_inlet), + doc="Levelized Cost of Treatment", + ) + + +def initialize_costing(m): + m.fs.costing.initialize() + + +def display_metrics_results(m): + print("----------Levelized costs----------") + LCOT = value( + pyunits.convert( + m.fs.costing.LCOT, to_units=m.fs.costing.base_currency / pyunits.m**3 + ) + ) + print(f"Levelized Cost of Treatment: {LCOT:.4f} $/m3 of feed") + LCOW = value( + pyunits.convert( + m.fs.costing.LCOW, to_units=m.fs.costing.base_currency / pyunits.m**3 + ) + ) + print(f"Levelized Cost of Water: {LCOW:.4f} $/m3 of product") + + print("----------Capital costs----------") + DCC_normalized = value( + pyunits.convert( + (m.fs.suboxicASM.costing.capital_cost) + / m.fs.costing.TIC + / m.fs.feed.properties[0].flow_vol, + to_units=m.fs.costing.base_currency / (pyunits.m**3 / pyunits.day), + ) + ) + print(f"Normalized direct capital costs: {DCC_normalized:.4f} $/(m3/day)") + ICC_normalized = value( + pyunits.convert( + m.fs.costing.total_capital_cost / m.fs.feed.properties[0].flow_vol, + to_units=m.fs.costing.base_currency / (pyunits.m**3 / pyunits.day), + ) + ) + print(f"Normalized total capital costs: {ICC_normalized:.4f} $/(m3/day)") + + print("----------Operating costs----------") + FMC_normalized = value( + pyunits.convert( + m.fs.costing.maintenance_cost / m.fs.costing.total_capital_cost, + to_units=1 / pyunits.a, + ) + ) + print(f"Normalized maintenance costs: {FMC_normalized:.4f} 1/year") + EC_normalized = value( + pyunits.convert( + m.fs.costing.aggregate_flow_costs["electricity"] + / m.fs.costing.annual_water_inlet, + to_units=m.fs.costing.base_currency / pyunits.m**3, + ) + ) + print(f"Normalized electricity cost: {EC_normalized:.4f} $/m3 of feed") + + print("----------Performance metrics----------") + volumetric_recovery = value( + m.fs.product_H2O.properties[0].flow_vol / m.fs.feed.properties[0].flow_vol + ) + print(f"Water recovery: {volumetric_recovery:.4f} m3 of product/m3 of feed") + BODR_normalized = value( + pyunits.convert( + 1 + - m.fs.product_H2O.properties[0].flow_mass_comp["bod"] + / m.fs.feed.properties[0].flow_mass_comp["bod"], + to_units=pyunits.dimensionless, + ) + ) + print(f"BOD removal: {BODR_normalized:.4f} dimensionless") + TSSR_normalized = value( + pyunits.convert( + 1 + - m.fs.product_H2O.properties[0].flow_mass_comp["tss"] + / m.fs.feed.properties[0].flow_mass_comp["tss"], + to_units=pyunits.dimensionless, + ) + ) + print(f"TSS removal: {TSSR_normalized:.4f} dimensionless") + TKNR_normalized = value( + pyunits.convert( + 1 + - (m.fs.product_H2O.properties[0].flow_mass_comp["tkn"]) + / (m.fs.feed.properties[0].flow_mass_comp["tkn"]), + to_units=pyunits.dimensionless, + ) + ) + print(f"TKN removal: {TKNR_normalized:.4f} dimensionless") + TP_normalized = value( + pyunits.convert( + 1 + - m.fs.product_H2O.properties[0].flow_mass_comp["phosphorus"] + / m.fs.feed.properties[0].flow_mass_comp["phosphorus"], + to_units=pyunits.dimensionless, + ) + ) + print(f"Total P removal: {TP_normalized:.4f} dimensionless") + + print("----------Energy intensity----------") + SEC = value( + pyunits.convert( + m.fs.costing.aggregate_flow_electricity / m.fs.feed.properties[0].flow_vol, + to_units=pyunits.kWh / pyunits.m**3, + ) + ) + print(f"Specific electricity consumption: {SEC:.4f} kWh/m3 of feed") + + +def display_additional_results(m): + print("----------Outlets----------") + product_H2O_flow = value( + pyunits.convert( + m.fs.product_H2O.properties[0].flow_vol, + to_units=pyunits.m**3 / pyunits.hr, + ) + ) + print(f"H2O outlet flow: {product_H2O_flow:.4f} m3/h") + product_H2O_BOD = value( + pyunits.convert( + m.fs.product_H2O.properties[0].conc_mass_comp["bod"], + to_units=pyunits.g / pyunits.L, + ) + ) + print(f"H2O outlet BOD conc: {product_H2O_BOD:.4f} g/L") + product_H2O_TSS = value( + pyunits.convert( + m.fs.product_H2O.properties[0].conc_mass_comp["tss"], + to_units=pyunits.g / pyunits.L, + ) + ) + print(f"H2O outlet TSS conc: {product_H2O_TSS:.4f} g/L") + product_H2O_TKN = value( + pyunits.convert( + m.fs.product_H2O.properties[0].conc_mass_comp["tkn"], + to_units=pyunits.g / pyunits.L, + ) + ) + print(f"H2O outlet TKN conc: {product_H2O_TKN:.4f} g/L") + product_H2O_TP = value( + pyunits.convert( + m.fs.product_H2O.properties[0].conc_mass_comp["phosphorus"], + to_units=pyunits.g / pyunits.L, + ) + ) + print(f"H2O outlet total P conc: {product_H2O_TP:.4f} g/L") + + print("----------Capital costs----------") + total_capital_costs = value(m.fs.costing.total_capital_cost) + print(f"Total capital costs: {total_capital_costs:.4f} $") + suboxicASM_capital_costs = value(m.fs.suboxicASM.costing.capital_cost) + print(f"Suboxic ASM capital costs: {suboxicASM_capital_costs:.4f} $") + + print("----------Operating costs----------") + total_operating_costs = value(m.fs.costing.total_operating_cost) / 1e6 + print(f"Total operating costs: {total_operating_costs:.4f} $M/year") + fixed_operating_costs = value(m.fs.costing.total_fixed_operating_cost) / 1e6 + print(f"Fixed operating costs: {fixed_operating_costs:.4f} $M/year") + electricity_operating_costs = ( + value(m.fs.costing.aggregate_flow_costs["electricity"]) / 1e3 + ) + print(f"Electricity operating costs: {electricity_operating_costs:.4f} $k/year") + + +if __name__ == "__main__": + m, results = main() diff --git a/watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml b/watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml new file mode 100644 index 0000000000..f616f97b2c --- /dev/null +++ b/watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml @@ -0,0 +1,46 @@ +base_currency: USD_2020 +base_period: year +defined_flows: + electricity: + value: 0.08 + units: USD_2020/kWh +global_parameters: + plant_lifetime: + value: 30 + units: year + utilization_factor: + value: 0.85 + units: dimensionless + land_cost_percent_FCI: + value: 0 + units: dimensionless + working_capital_percent_FCI: + value: 0 + units: dimensionless + salaries_percent_FCI: + value: 0 + units: 1/year + benefit_percent_of_salary: + value: 0 + units: dimensionless + maintenance_costs_percent_FCI: + value: 0.02 + units: 1/year + laboratory_fees_percent_FCI: + value: 0 + units: 1/year + insurance_and_taxes_percent_FCI: + value: 0 + units: 1/year + wacc: + # Weighted Average Cost of Capital (WACC) + value: 0.05 + units: dimensionless + TPEC: + # Total Purchased Equipment Cost + value: 3.4 + units: dimensionless + TIC: + # Total Installed Cost + value: 2 + units: dimensionless \ No newline at end of file From 258f6fe4d56bc0a718ae08ee4c21833a7dd437f9 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Tue, 11 Apr 2023 15:07:10 -0400 Subject: [PATCH 02/21] delete redundant files --- .../flowsheets/storage_tank/storage_tank.py | 340 ------------------ ...boxic_activated_sludge_process_global.yaml | 46 --- 2 files changed, 386 deletions(-) delete mode 100644 watertap/examples/flowsheets/storage_tank/storage_tank.py delete mode 100644 watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml diff --git a/watertap/examples/flowsheets/storage_tank/storage_tank.py b/watertap/examples/flowsheets/storage_tank/storage_tank.py deleted file mode 100644 index 9ae3a55eb4..0000000000 --- a/watertap/examples/flowsheets/storage_tank/storage_tank.py +++ /dev/null @@ -1,340 +0,0 @@ -################################################################################# -# WaterTAP Copyright (c) 2020-2023, 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 os -import idaes.logger as idaeslog -from pyomo.environ import ( - ConcreteModel, - Expression, - value, - TransformationFactory, - units as pyunits, -) -from pyomo.network import Arc, SequentialDecomposition -from pyomo.util.check_units import assert_units_consistent - -from idaes.core import FlowsheetBlock -from idaes.core.solvers import get_solver -import idaes.core.util.scaling as iscale -from idaes.models.unit_models import Product -from idaes.core import UnitModelCostingBlock - -from watertap.core.wt_database import Database -import watertap.core.zero_order_properties as prop_ZO -from watertap.core.util.initialization import assert_degrees_of_freedom, check_solve -from watertap.unit_models.zero_order import ( - FeedZO, - StorageTankZO, -) -from watertap.core.zero_order_costing import ZeroOrderCosting - -# Set up logger -_log = idaeslog.getLogger(__name__) - - -def main(): - m = build() - - set_operating_conditions(m) - assert_degrees_of_freedom(m, 0) - assert_units_consistent(m) - - initialize_system(m) - - results = solve(m, checkpoint="solve flowsheet after initializing system") - display_results(m) - - add_costing(m) - initialize_costing(m) - assert_degrees_of_freedom(m, 0) - assert_units_consistent(m) - - results = solve(m, checkpoint="solve flowsheet after costing") - - display_metrics_results(m) - display_additional_results(m) - return m, results - - -def build(): - m = ConcreteModel() - m.db = Database() - - m.fs = FlowsheetBlock(dynamic=False) - m.fs.prop = prop_ZO.WaterParameterBlock( - solute_list=["bod", "tss", "tkn", "phosphorus"] - ) - - # components - m.fs.feed = FeedZO(property_package=m.fs.prop) - m.fs.product_H2O = Product(property_package=m.fs.prop) - - m.fs.suboxicASM = StorageTankZO(property_package=m.fs.prop, database=m.db) - - # connections - m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.suboxicASM.inlet) - m.fs.s02 = Arc(source=m.fs.suboxicASM.outlet, destination=m.fs.product_H2O.inlet) - - TransformationFactory("network.expand_arcs").apply_to(m) - - # scale - iscale.calculate_scaling_factors(m) - return m - - -def set_operating_conditions(m): - # feed - flow_vol = 10 * pyunits.Mgallons / pyunits.day - conc_bod = 314 * pyunits.mg / pyunits.L - conc_tss = 336 * pyunits.mg / pyunits.L - conc_tkn = 52 * pyunits.mg / pyunits.L - conc_phosphorus = 6 * pyunits.mg / pyunits.L - - m.fs.feed.flow_vol[0].fix(flow_vol) - m.fs.feed.conc_mass_comp[0, "bod"].fix(conc_bod) - m.fs.feed.conc_mass_comp[0, "tss"].fix(conc_tss) - m.fs.feed.conc_mass_comp[0, "tkn"].fix(conc_tkn) - m.fs.feed.conc_mass_comp[0, "phosphorus"].fix(conc_phosphorus) - solve(m.fs.feed, checkpoint="solve feed block") - - # suboxicASM - m.fs.suboxicASM.load_parameters_from_database(use_default_removal=True) - - -def initialize_system(m): - seq = SequentialDecomposition() - seq.options.tear_set = [] - seq.options.iterLim = 1 - seq.run(m, lambda u: u.initialize()) - - -def solve(blk, solver=None, checkpoint=None, tee=False, fail_flag=True): - if solver is None: - solver = get_solver() - results = solver.solve(blk, tee=tee) - check_solve(results, checkpoint=checkpoint, logger=_log, fail_flag=fail_flag) - return results - - -def display_results(m): - unit_list = ["feed", "suboxicASM"] - for u in unit_list: - m.fs.component(u).report() - - -def add_costing(m): - source_file = os.path.join( - os.path.dirname(os.path.abspath(__file__)), - "suboxic_activated_sludge_process_global.yaml", - ) - m.fs.costing = ZeroOrderCosting(case_study_definition=source_file) - # typing aid - costing_kwargs = {"flowsheet_costing_block": m.fs.costing} - m.fs.suboxicASM.costing = UnitModelCostingBlock(**costing_kwargs) - - m.fs.costing.cost_process() - m.fs.costing.add_electricity_intensity(m.fs.product_H2O.properties[0].flow_vol) - m.fs.costing.add_LCOW(m.fs.product_H2O.properties[0].flow_vol) - - # other levelized costs - m.fs.costing.annual_water_inlet = Expression( - expr=m.fs.costing.utilization_factor - * pyunits.convert( - m.fs.feed.properties[0].flow_vol, - to_units=pyunits.m**3 / m.fs.costing.base_period, - ) - ) - - m.fs.costing.annual_water_production = Expression( - expr=m.fs.costing.utilization_factor - * pyunits.convert( - m.fs.product_H2O.properties[0].flow_vol, - to_units=pyunits.m**3 / m.fs.costing.base_period, - ) - ) - - m.fs.costing.total_annualized_cost = Expression( - expr=( - m.fs.costing.total_capital_cost * m.fs.costing.capital_recovery_factor - + m.fs.costing.total_operating_cost - ) - ) - - m.fs.costing.LCOT = Expression( - expr=(m.fs.costing.total_annualized_cost / m.fs.costing.annual_water_inlet), - doc="Levelized Cost of Treatment", - ) - - -def initialize_costing(m): - m.fs.costing.initialize() - - -def display_metrics_results(m): - print("----------Levelized costs----------") - LCOT = value( - pyunits.convert( - m.fs.costing.LCOT, to_units=m.fs.costing.base_currency / pyunits.m**3 - ) - ) - print(f"Levelized Cost of Treatment: {LCOT:.4f} $/m3 of feed") - LCOW = value( - pyunits.convert( - m.fs.costing.LCOW, to_units=m.fs.costing.base_currency / pyunits.m**3 - ) - ) - print(f"Levelized Cost of Water: {LCOW:.4f} $/m3 of product") - - print("----------Capital costs----------") - DCC_normalized = value( - pyunits.convert( - (m.fs.suboxicASM.costing.capital_cost) - / m.fs.costing.TIC - / m.fs.feed.properties[0].flow_vol, - to_units=m.fs.costing.base_currency / (pyunits.m**3 / pyunits.day), - ) - ) - print(f"Normalized direct capital costs: {DCC_normalized:.4f} $/(m3/day)") - ICC_normalized = value( - pyunits.convert( - m.fs.costing.total_capital_cost / m.fs.feed.properties[0].flow_vol, - to_units=m.fs.costing.base_currency / (pyunits.m**3 / pyunits.day), - ) - ) - print(f"Normalized total capital costs: {ICC_normalized:.4f} $/(m3/day)") - - print("----------Operating costs----------") - FMC_normalized = value( - pyunits.convert( - m.fs.costing.maintenance_cost / m.fs.costing.total_capital_cost, - to_units=1 / pyunits.a, - ) - ) - print(f"Normalized maintenance costs: {FMC_normalized:.4f} 1/year") - EC_normalized = value( - pyunits.convert( - m.fs.costing.aggregate_flow_costs["electricity"] - / m.fs.costing.annual_water_inlet, - to_units=m.fs.costing.base_currency / pyunits.m**3, - ) - ) - print(f"Normalized electricity cost: {EC_normalized:.4f} $/m3 of feed") - - print("----------Performance metrics----------") - volumetric_recovery = value( - m.fs.product_H2O.properties[0].flow_vol / m.fs.feed.properties[0].flow_vol - ) - print(f"Water recovery: {volumetric_recovery:.4f} m3 of product/m3 of feed") - BODR_normalized = value( - pyunits.convert( - 1 - - m.fs.product_H2O.properties[0].flow_mass_comp["bod"] - / m.fs.feed.properties[0].flow_mass_comp["bod"], - to_units=pyunits.dimensionless, - ) - ) - print(f"BOD removal: {BODR_normalized:.4f} dimensionless") - TSSR_normalized = value( - pyunits.convert( - 1 - - m.fs.product_H2O.properties[0].flow_mass_comp["tss"] - / m.fs.feed.properties[0].flow_mass_comp["tss"], - to_units=pyunits.dimensionless, - ) - ) - print(f"TSS removal: {TSSR_normalized:.4f} dimensionless") - TKNR_normalized = value( - pyunits.convert( - 1 - - (m.fs.product_H2O.properties[0].flow_mass_comp["tkn"]) - / (m.fs.feed.properties[0].flow_mass_comp["tkn"]), - to_units=pyunits.dimensionless, - ) - ) - print(f"TKN removal: {TKNR_normalized:.4f} dimensionless") - TP_normalized = value( - pyunits.convert( - 1 - - m.fs.product_H2O.properties[0].flow_mass_comp["phosphorus"] - / m.fs.feed.properties[0].flow_mass_comp["phosphorus"], - to_units=pyunits.dimensionless, - ) - ) - print(f"Total P removal: {TP_normalized:.4f} dimensionless") - - print("----------Energy intensity----------") - SEC = value( - pyunits.convert( - m.fs.costing.aggregate_flow_electricity / m.fs.feed.properties[0].flow_vol, - to_units=pyunits.kWh / pyunits.m**3, - ) - ) - print(f"Specific electricity consumption: {SEC:.4f} kWh/m3 of feed") - - -def display_additional_results(m): - print("----------Outlets----------") - product_H2O_flow = value( - pyunits.convert( - m.fs.product_H2O.properties[0].flow_vol, - to_units=pyunits.m**3 / pyunits.hr, - ) - ) - print(f"H2O outlet flow: {product_H2O_flow:.4f} m3/h") - product_H2O_BOD = value( - pyunits.convert( - m.fs.product_H2O.properties[0].conc_mass_comp["bod"], - to_units=pyunits.g / pyunits.L, - ) - ) - print(f"H2O outlet BOD conc: {product_H2O_BOD:.4f} g/L") - product_H2O_TSS = value( - pyunits.convert( - m.fs.product_H2O.properties[0].conc_mass_comp["tss"], - to_units=pyunits.g / pyunits.L, - ) - ) - print(f"H2O outlet TSS conc: {product_H2O_TSS:.4f} g/L") - product_H2O_TKN = value( - pyunits.convert( - m.fs.product_H2O.properties[0].conc_mass_comp["tkn"], - to_units=pyunits.g / pyunits.L, - ) - ) - print(f"H2O outlet TKN conc: {product_H2O_TKN:.4f} g/L") - product_H2O_TP = value( - pyunits.convert( - m.fs.product_H2O.properties[0].conc_mass_comp["phosphorus"], - to_units=pyunits.g / pyunits.L, - ) - ) - print(f"H2O outlet total P conc: {product_H2O_TP:.4f} g/L") - - print("----------Capital costs----------") - total_capital_costs = value(m.fs.costing.total_capital_cost) - print(f"Total capital costs: {total_capital_costs:.4f} $") - suboxicASM_capital_costs = value(m.fs.suboxicASM.costing.capital_cost) - print(f"Suboxic ASM capital costs: {suboxicASM_capital_costs:.4f} $") - - print("----------Operating costs----------") - total_operating_costs = value(m.fs.costing.total_operating_cost) / 1e6 - print(f"Total operating costs: {total_operating_costs:.4f} $M/year") - fixed_operating_costs = value(m.fs.costing.total_fixed_operating_cost) / 1e6 - print(f"Fixed operating costs: {fixed_operating_costs:.4f} $M/year") - electricity_operating_costs = ( - value(m.fs.costing.aggregate_flow_costs["electricity"]) / 1e3 - ) - print(f"Electricity operating costs: {electricity_operating_costs:.4f} $k/year") - - -if __name__ == "__main__": - m, results = main() diff --git a/watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml b/watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml deleted file mode 100644 index f616f97b2c..0000000000 --- a/watertap/examples/flowsheets/storage_tank/suboxic_activated_sludge_process_global.yaml +++ /dev/null @@ -1,46 +0,0 @@ -base_currency: USD_2020 -base_period: year -defined_flows: - electricity: - value: 0.08 - units: USD_2020/kWh -global_parameters: - plant_lifetime: - value: 30 - units: year - utilization_factor: - value: 0.85 - units: dimensionless - land_cost_percent_FCI: - value: 0 - units: dimensionless - working_capital_percent_FCI: - value: 0 - units: dimensionless - salaries_percent_FCI: - value: 0 - units: 1/year - benefit_percent_of_salary: - value: 0 - units: dimensionless - maintenance_costs_percent_FCI: - value: 0.02 - units: 1/year - laboratory_fees_percent_FCI: - value: 0 - units: 1/year - insurance_and_taxes_percent_FCI: - value: 0 - units: 1/year - wacc: - # Weighted Average Cost of Capital (WACC) - value: 0.05 - units: dimensionless - TPEC: - # Total Purchased Equipment Cost - value: 3.4 - units: dimensionless - TIC: - # Total Installed Cost - value: 2 - units: dimensionless \ No newline at end of file From d7c1d78ad68f82f48f3637ea11c7b7d7ebc65727 Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Mon, 9 Jun 2025 18:37:12 -0400 Subject: [PATCH 03/21] first draft of alternative activated sludge process (reactors only) --- .../activated_sludge_reactor_train.py | 393 ++++++++++++++++++ 1 file changed, 393 insertions(+) create mode 100644 watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py diff --git a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py b/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py new file mode 100644 index 0000000000..6dc789b886 --- /dev/null +++ b/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py @@ -0,0 +1,393 @@ +################################################################################# +# 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/" +################################################################################# +""" +Example of activated sludge process model. +Initial layout: + * 5 reactors + * R1, R3, R5: anoxic + * R2 and R4: aerobic + * 2 Mixers: + * M1 mixes feed and R5 split fraction, outlet to R1 inlet + * M2 mixes R1 outlet, R2 outlet, outlet to R3 inlet + * 1 Splitter: + * separates R5 outlet into effluent, M1 inlet, and R2 inlet + + +Unit operations are modeled as follows: + + * Anoxic reactors as standard CSTRs (CSTR) + * Aerobic reactors as AerationTanks + +Note that a pressure changer is required in the recycle stream to ensure the +pressure inside the recycle loop is bounded. As the inlet Mixer uses a pressure +minimization constraint and there is no pressure drop in the reactors, if pressure +is not specified at some point within the recycle loop then it becomes unbounded. + +""" + +__author__ = "Adam Atia" + +import pyomo.environ as pyo +from pyomo.network import Arc, SequentialDecomposition + +from idaes.core import FlowsheetBlock +from idaes.models.unit_models import ( + CSTR, + Feed, + Mixer, + Separator, + PressureChanger, + Product, + MomentumMixingType +) +from idaes.models.unit_models.separator import SplittingType +from watertap.core.solvers import get_solver +from idaes.core.util.model_statistics import degrees_of_freedom +import idaes.logger as idaeslog +import idaes.core.util.scaling as iscale +from idaes.core.util.tables import ( + create_stream_table_dataframe, + stream_table_dataframe_to_string, +) +from idaes.core.util.initialization import propagate_state + +from watertap.unit_models import AerationTank + +from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( + ASM1ParameterBlock, +) +from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( + ASM1ReactionParameterBlock, +) + +from watertap.core.util.initialization import check_solve, interval_initializer + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +def build_flowsheet(): + m = pyo.ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ASM1ParameterBlock() + m.fs.rxn_props = ASM1ReactionParameterBlock(property_package=m.fs.props) + # Feed water stream + m.fs.feed = Feed(property_package=m.fs.props) + + # Mixer for feed water and recycled sludge + m.fs.M1 = Mixer(property_package=m.fs.props, + inlet_list=["feed", "recycle"], + # momentum_mixing_type=MomentumMixingType.equality + ) + # m.fs.M1.pressure_equality_constraints[0,2].deactivate() + # First reactor (anoxic) + m.fs.R1 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) + + # Second reactor (aerobic) + m.fs.R2 = AerationTank(property_package=m.fs.props, reaction_package=m.fs.rxn_props) + + + # Mixer for R1 and R2 outlets + m.fs.M2= Mixer(property_package=m.fs.props, inlet_list=["R1_outlet", "R2_outlet"]) + + + # Third reactor (anoxic) + m.fs.R3 = CSTR( + property_package=m.fs.props, reaction_package=m.fs.rxn_props + ) + # Fourth reactor (aerobic) - CSTR with injection + m.fs.R4 = AerationTank( + property_package=m.fs.props, reaction_package=m.fs.rxn_props + ) + # Fifth reactor (anoxic) + m.fs.R5 = CSTR( + property_package=m.fs.props, reaction_package=m.fs.rxn_props + ) + m.fs.S1 = Separator( + property_package=m.fs.props, outlet_list=["effluent","M1_inlet", "R2_inlet"] + ) + # Clarifier + # # TODO: Replace with more detailed model when available + # m.fs.CL1 = Separator( + # property_package=m.fs.props, + # outlet_list=["underflow", "effluent"], + # split_basis=SplittingType.componentFlow, + # ) + # # Sludge purge splitter + # m.fs.SP6 = Separator( + # property_package=m.fs.props, + # outlet_list=["recycle", "waste"], + # split_basis=SplittingType.totalFlow, + # ) + # # Mixing sludge recycle and R5 underflow + # m.fs.MX6 = Mixer(property_package=m.fs.props, inlet_list=["clarifier", "reactor"]) + + # Product Blocks + m.fs.Treated = Product(property_package=m.fs.props) + # Recycle pressure changer - use a simple isothermal unit for now + # m.fs.P1 = PressureChanger(property_package=m.fs.props) + + # Link units + m.fs.feed_to_m1 = Arc(source=m.fs.feed.outlet, destination=m.fs.M1.feed) + m.fs.m1_to_r1 = Arc(source=m.fs.M1.outlet, destination=m.fs.R1.inlet) + m.fs.r1_to_m2 = Arc(source=m.fs.R1.outlet, destination=m.fs.M2.R1_outlet) + m.fs.r2_to_m2 = Arc(source=m.fs.R2.outlet, destination=m.fs.M2.R2_outlet) + m.fs.m2_to_r3 = Arc(source=m.fs.M2.outlet, destination=m.fs.R3.inlet) + m.fs.r3_to_r4 = Arc(source=m.fs.R3.outlet, destination=m.fs.R4.inlet) + m.fs.r4_to_r5 = Arc(source=m.fs.R4.outlet, destination=m.fs.R5.inlet) + m.fs.r5_to_s1 = Arc(source=m.fs.R5.outlet, destination=m.fs.S1.inlet) + m.fs.s1_to_effluent = Arc(source=m.fs.S1.effluent, destination=m.fs.Treated.inlet) + m.fs.s1_to_m1 = Arc(source=m.fs.S1.M1_inlet, destination=m.fs.M1.recycle) + m.fs.s1_to_r2 = Arc(source=m.fs.S1.R2_inlet, destination=m.fs.R2.inlet) + # m.fs.stream14 = Arc(source=m.fs.MX6.outlet, destination=m.fs.P1.inlet) + # m.fs.stream15 = Arc(source=m.fs.P1.outlet, destination=m.fs.MX1.recycle) + pyo.TransformationFactory("network.expand_arcs").apply_to(m) + + return m + +def set_operating_conditions(m): + # Feed Water Conditions + m.fs.feed.flow_vol.fix(18446 * pyo.units.m**3 / pyo.units.day) + m.fs.feed.temperature.fix(298.15 * pyo.units.K) + m.fs.feed.pressure.fix(1 * pyo.units.atm) + m.fs.feed.conc_mass_comp[0, "S_I"].fix(30 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "S_S"].fix(69.5 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_I"].fix(51.2 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_S"].fix(202.32 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_BH"].fix(28.17 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_BA"].fix(0 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_P"].fix(0 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "S_O"].fix(0 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "S_NO"].fix(0 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "S_NH"].fix(31.56 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "S_ND"].fix(6.95 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_ND"].fix(10.59 * pyo.units.g / pyo.units.m**3) + m.fs.feed.alkalinity.fix(7 * pyo.units.mol / pyo.units.m**3) + + # Reactor sizing + m.fs.R1.volume.fix(1000 * pyo.units.m**3) + m.fs.R2.volume.fix(1000 * pyo.units.m**3) + m.fs.R3.volume.fix(1333 * pyo.units.m**3) + m.fs.R4.volume.fix(1333 * pyo.units.m**3) + m.fs.R5.volume.fix(1333 * pyo.units.m**3) + + # Injection rates to Reactors 2 and 4 + for j in m.fs.props.component_list: + if j != "S_O": + # All components except S_O have no injection + m.fs.R2.injection[:, :, j].fix(0) + m.fs.R4.injection[:, :, j].fix(0) + # Then set injections rates for O2 + m.fs.R2.outlet.conc_mass_comp[:, "S_O"].fix(1.72e-3) + m.fs.R4.outlet.conc_mass_comp[:, "S_O"].fix(2.43e-3) + + # Set fraction of outflow from reactor 5 that goes to recycle + m.fs.S1.split_fraction[:, "M1_inlet"].fix(0.5) + + # Set fraction of outflow from reactor 5 that goes to effluent + m.fs.S1.split_fraction[:, "effluent"].fix(0.4) + + # # Outlet pressure from recycle pump + # m.fs.P1.outlet.pressure.fix(101325) + + # Check degrees of freedom + print("DOF = ", degrees_of_freedom(m)) + assert degrees_of_freedom(m) == 0 + +def scale_flowsheet(m): + # Apply scaling + for var in m.fs.component_data_objects(pyo.Var, descend_into=True): + if "flow_vol" in var.name: + iscale.set_scaling_factor(var, 1e1) + if "temperature" in var.name: + iscale.set_scaling_factor(var, 1e-1) + if "pressure" in var.name: + iscale.set_scaling_factor(var, 1e-6) + if "conc_mass_comp" in var.name: + iscale.set_scaling_factor(var, 1e2) + iscale.calculate_scaling_factors(m.fs) + +def init_and_propagate(blk,arc=None,source=None,destination=None): + blk.initialize() + if arc is not None: + propagate_state(arc) + elif source is not None: + propagate_state(source=source,destination=destination) + + +def initialize_flowsheet(m): + # Initialize flowsheet + # interval_initializer(m) + m.fs.feed.initialize() + propagate_state(m.fs.feed_to_m1) + propagate_state(source=m.fs.feed.outlet, destination=m.fs.M1.recycle) + + m.fs.M1.initialize() + propagate_state(m.fs.m1_to_r1) + + m.fs.R1.initialize() + propagate_state(m.fs.r1_to_m2) + propagate_state(source=m.fs.R1.outlet, destination=m.fs.M2.R2_outlet) + + m.fs.M2.initialize() + propagate_state(m.fs.m2_to_r3) + + m.fs.R3.initialize() + propagate_state(m.fs.r3_to_r4) + + m.fs.R4.initialize() + propagate_state(m.fs.r4_to_r5) + + m.fs.R5.initialize() + propagate_state(m.fs.r5_to_s1) + + m.fs.S1.initialize() + propagate_state(m.fs.s1_to_effluent) + propagate_state(m.fs.s1_to_m1) + propagate_state(m.fs.s1_to_r2) + + m.fs.R2.initialize() + propagate_state(m.fs.r2_to_m2) + + m.fs.M1.initialize() + propagate_state(m.fs.m1_to_r1) + + m.fs.R1.initialize() + propagate_state(m.fs.r1_to_m2) + + m.fs.M2.initialize() + propagate_state(m.fs.m2_to_r3) + + m.fs.R3.initialize() + propagate_state(m.fs.r3_to_r4) + + m.fs.R4.initialize() + propagate_state(m.fs.r4_to_r5) + + m.fs.R5.initialize() + propagate_state(m.fs.r5_to_s1) + + m.fs.S1.initialize() + propagate_state(m.fs.s1_to_effluent) + propagate_state(m.fs.s1_to_m1) + propagate_state(m.fs.s1_to_r2) + + m.fs.Treated.initialize() + + interval_initializer(m) + + + # Apply sequential decomposition - 1 iteration should suffice + seq = SequentialDecomposition() + seq.options.select_tear_method = "heuristic" + seq.options.tear_method = "Wegstein" + seq.options.iterLim = 1 + + 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") + # order = seq.calculation_order(G) + # for o in heuristic_tear_set: + # print(o.name) + # for o in order: + # print(o[0].name) + + # # Initial guesses for flow into first reactor + # tear_guesses = { + # "flow_vol": {0: 1.0675}, + # "conc_mass_comp": { + # (0, "S_I"): 0.03, + # (0, "S_S"): 0.0146, + # (0, "X_I"): 1.149, + # (0, "X_S"): 0.0893, + # (0, "X_BH"): 2.542, + # (0, "X_BA"): 0.149, + # (0, "X_P"): 0.448, + # (0, "S_O"): 3.93e-4, + # (0, "S_NO"): 8.32e-3, + # (0, "S_NH"): 7.7e-3, + # (0, "S_ND"): 1.94e-3, + # (0, "X_ND"): 5.62e-3, + # }, + # "alkalinity": {0: 5e-3}, + # "temperature": {0: 298.15}, + # "pressure": {0: 101325}, + # } + + # # Pass the tear_guess to the SD tool + # seq.set_guesses_for(m.fs.R1.inlet, tear_guesses) + + def function(unit): + unit.initialize(outlvl=idaeslog.DEBUG) + + seq.run(m, function) + +def solve_flowsheet(m): + # Solve overall flowsheet to close recycle loop + solver = get_solver() + results = solver.solve(m) + check_solve(results, checkpoint="closing recycle", logger=_log, fail_flag=True) + + # Switch to fixed KLa in R3 and R4 (S_O concentration is controlled in R5) + m.fs.R2.KLa.fix(10) + m.fs.R4.KLa.fix(10) + m.fs.R2.outlet.conc_mass_comp[:, "S_O"].unfix() + m.fs.R4.outlet.conc_mass_comp[:, "S_O"].unfix() + + # Resolve with controls in place + results = solver.solve(m, tee=True) + check_solve( + results, + checkpoint="re-solve with controls in place", + logger=_log, + fail_flag=True, + ) + + return m, results + + +if __name__ == "__main__": + # This method builds and runs a steady state activated sludge + # flowsheet. + # m, results = build_flowsheet() + m = build_flowsheet() + set_operating_conditions(m) + scale_flowsheet(m) + + initialize_flowsheet(m) + + scale_flowsheet(m) + + m, res = solve_flowsheet(m) + + + stream_table = create_stream_table_dataframe( + { + "Feed": m.fs.feed.outlet, + "M1": m.fs.M1.outlet, + "R1": m.fs.R1.outlet, + "M2": m.fs.M2.outlet, + "R2": m.fs.R2.outlet, + "R3": m.fs.R3.outlet, + "R4": m.fs.R4.outlet, + "R5": m.fs.R5.outlet, + "S1 to M1": m.fs.S1.M1_inlet, + "S1 to R2": m.fs.S1.R2_inlet, + "Effluent": m.fs.Treated.inlet + + }, + time_point=0, + ) + print(stream_table_dataframe_to_string(stream_table)) From 0752d35bca2e9e13c02d3b8b913c28ff57b1cf9a Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Mon, 9 Jun 2025 18:48:33 -0400 Subject: [PATCH 04/21] make COD_to_SS a param on ASM1 for consistency --- .../unit_specific/activated_sludge/asm1_properties.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py index 99958101a2..d91dff088f 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py @@ -144,7 +144,7 @@ def build(self): domain=pyo.PositiveReals, doc="Mass fraction of N per COD in particulates, i_xp", ) - self.COD_to_SS = pyo.Var( + self.COD_to_SS = pyo.Param( initialize=0.75, units=pyo.units.dimensionless, domain=pyo.PositiveReals, From 8bd13db905c4d2524c638fb8701f4fca77cabc19 Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Mon, 9 Jun 2025 19:06:22 -0400 Subject: [PATCH 05/21] tinkering with scaling --- .../activated_sludge_reactor_train.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py b/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py index 6dc789b886..c5be8b4032 100644 --- a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py +++ b/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py @@ -88,8 +88,9 @@ def build_flowsheet(): # Mixer for feed water and recycled sludge m.fs.M1 = Mixer(property_package=m.fs.props, inlet_list=["feed", "recycle"], - # momentum_mixing_type=MomentumMixingType.equality + momentum_mixing_type=MomentumMixingType.none ) + m.fs.M1.outlet.pressure.fix() # m.fs.M1.pressure_equality_constraints[0,2].deactivate() # First reactor (anoxic) m.fs.R1 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) @@ -99,8 +100,10 @@ def build_flowsheet(): # Mixer for R1 and R2 outlets - m.fs.M2= Mixer(property_package=m.fs.props, inlet_list=["R1_outlet", "R2_outlet"]) - + m.fs.M2= Mixer(property_package=m.fs.props, inlet_list=["R1_outlet", "R2_outlet"], + momentum_mixing_type=MomentumMixingType.none + ) + m.fs.M2.outlet.pressure.fix() # Third reactor (anoxic) m.fs.R3 = CSTR( @@ -213,7 +216,7 @@ def scale_flowsheet(m): if "temperature" in var.name: iscale.set_scaling_factor(var, 1e-1) if "pressure" in var.name: - iscale.set_scaling_factor(var, 1e-6) + iscale.set_scaling_factor(var, 1e-4) if "conc_mass_comp" in var.name: iscale.set_scaling_factor(var, 1e2) iscale.calculate_scaling_factors(m.fs) From f9048ec7e95b908458ab144677ab554d36864ae0 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Tue, 10 Jun 2025 16:13:08 -0400 Subject: [PATCH 06/21] ASM3 development --- .../activated_sludge/asm3_properties.py | 643 ++++++++++++++++++ .../activated_sludge/asm3_reactions.py | 603 ++++++++++++++++ 2 files changed, 1246 insertions(+) create mode 100644 watertap/property_models/unit_specific/activated_sludge/asm3_properties.py create mode 100644 watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py new file mode 100644 index 0000000000..5b403f4ea2 --- /dev/null +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py @@ -0,0 +1,643 @@ +################################################################################# +# 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/" +################################################################################# +""" +Thermophysical property package to be used in conjunction with ASM3 reactions. +""" + +# Import Pyomo libraries +import pyomo.environ as pyo + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + MaterialFlowBasis, + StateBlockData, + StateBlock, + MaterialBalanceType, + EnergyBalanceType, + LiquidPhase, + Component, + Solute, + Solvent, +) +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.initialization import fix_state_vars, revert_state_vars +from idaes.core.scaling import CustomScalerBase +import idaes.logger as idaeslog +from idaes.core.base.property_base import PhysicalParameterBlock +import idaes.core.util.scaling as iscale + +# Some more information about this module +__author__ = "Chenyu Wang" + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("ASM3ParameterBlock") +class ASM3ParameterData(PhysicalParameterBlock): + """ + Property Parameter Block Class + """ + + def build(self): + """ + Callable method for Block construction. + """ + super().build() + + self._state_block_class = ASM3StateBlock + + # Add Phase objects + self.Liq = LiquidPhase() + + # Add Component objects + self.H2O = Solvent() + + # Soluble Components + self.S_O = Solute(doc="Dissolved Oxygen, S_O") + self.S_I = Solute(doc="Inert soluble organic material, S_I") + self.S_S = Solute(doc="Readily biodegradable substrates, S_S") + self.S_NH = Solute( + doc="Ammonium plus ammonia nitrogen (NH4+ + NH3 nitrogen), S_NH" + ) + self.S_N2 = Solute(doc="Dinitrogen (N2), S_N2") + self.S_NOX = Solute( + doc="Nitrate plus nitrite nitrogen (NO3- + NO2- nitrogen), S_NOX" + ) + self.S_ALK = Component(doc="Alkalinity (HCO3-), S_ALK") + + # Particulate Components + self.X_I = Solute(doc="Inert particulate organic matter, X_I") + self.X_S = Solute(doc="Slowly biodegradable substrates, X_S") + self.X_H = Solute(doc="Heterotrophic organisms, X_H") + self.X_STO = Solute( + doc="A cell internal storage product of heterotrophic organisms, X_STO" + ) + self.X_A = Solute(doc="Nitrifying organisms, X_A") + self.X_TSS = Solute(doc="Total suspended solids, X_TSS") + + # Create sets for use across ASM models and associated unit models (e.g., thickener, dewaterer) + self.non_particulate_component_set = pyo.Set( + initialize=[ + "S_O", + "S_I", + "S_S", + "S_NH", + "S_N2", + "S_NOX", + "S_ALK", + "H2O", + ] + ) + self.particulate_component_set = pyo.Set( + initialize=["X_I", "X_S", "X_H", "X_STO", "X_A", "X_TSS"] + ) + + # Heat capacity of water + self.cp_mass = pyo.Param( + initialize=4182, + doc="Specific heat capacity of water", + units=pyo.units.J / pyo.units.kg / pyo.units.K, + ) + # Density of water + self.dens_mass = pyo.Param( + initialize=997, + doc="Density of water", + units=pyo.units.kg / pyo.units.m**3, + ) + + # Thermodynamic reference state + self.pressure_ref = pyo.Param( + within=pyo.PositiveReals, + mutable=True, + default=101325.0, + doc="Reference pressure", + units=pyo.units.Pa, + ) + self.temperature_ref = pyo.Param( + within=pyo.PositiveReals, + mutable=True, + default=298.15, + doc="Reference temperature", + units=pyo.units.K, + ) + + # Typical stoichiometric and composition parametersfor ASM3 + self.f_SI = pyo.Var( + initialize=0, + units=pyo.units.dimensionless, + domain=pyo.NonNegativeReals, + doc="Production of S_I in hydrolysis (g-COD-S_I / g-COD-X_S)", + ) + self.i_NSI = pyo.Var( + initialize=0.01, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="N content of S_I (g-N / g-COD-S_I)", + ) + self.i_NSS = pyo.Var( + initialize=0.03, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="N content of S_S (g-N / g-COD-S_S)", + ) + self.i_NXI = pyo.Var( + initialize=0.02, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="N content of X_I (g-N / g-COD-X_I)", + ) + self.i_NXS = pyo.Var( + initialize=0.04, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="N content of X_S (g-N / g-COD-X_S)", + ) + self.i_NBM = pyo.Var( + initialize=0.07, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="N content of biomass, X_H, X_A (g-N / g-COD-X_H or X_A)", + ) + self.i_TSXI = pyo.Var( + initialize=0.75, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="TSS to COD ratio for X_I (g-TSS / g-COD-X_I)", + ) + self.i_TSXS = pyo.Var( + initialize=0.75, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="TSS to COD ratio for X_S (g-TSS / g-COD-X_S)", + ) + self.i_TSBM = pyo.Var( + initialize=0.90, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="TSS to COD ratio for for biomass,X_H, X_A (g-TSS / g-COD-X_H or X_A)", + ) + self.i_TSSTO = pyo.Var( + initialize=0.60, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="TSS to COD ratio for X_STO based on PHB (g-TSS / g-X_STO)", + ) + # self.f_p = pyo.Var( + # initialize=0.08, + # units=pyo.units.dimensionless, + # domain=pyo.PositiveReals, + # doc="Fraction of biomass yielding particulate products, f_p", + # ) + # self.i_xb = pyo.Var( + # initialize=0.08, + # units=pyo.units.dimensionless, + # domain=pyo.PositiveReals, + # doc="Mass fraction of N per COD in biomass, i_xb", + # ) + # self.i_xp = pyo.Var( + # initialize=0.06, + # units=pyo.units.dimensionless, + # domain=pyo.PositiveReals, + # doc="Mass fraction of N per COD in particulates, i_xp", + # ) + # self.COD_to_SS = pyo.Var( + # initialize=0.75, + # units=pyo.units.dimensionless, + # domain=pyo.PositiveReals, + # doc="Conversion factor applied for TSS calculation", + # ) + # self.BOD5_factor = pyo.Param( + # ["raw", "effluent"], + # initialize={"raw": 0.65, "effluent": 0.25}, + # units=pyo.units.dimensionless, + # domain=pyo.PositiveReals, + # doc="Conversion factor for BOD5", + # ) + # Fix Vars that are treated as Params + for v in self.component_objects(pyo.Var): + v.fix() + + @classmethod + def define_metadata(cls, obj): + obj.add_properties( + { + "flow_vol": {"method": None}, + "pressure": {"method": None}, + "temperature": {"method": None}, + "conc_mass_comp": {"method": None}, + } + ) + # obj.define_custom_properties( + # { + # "alkalinity": {"method": None}, + # "TSS": {"method": "_TSS"}, + # "BOD5": {"method": "_BOD5"}, + # "TKN": {"method": "_TKN"}, + # "Total_N": {"method": "_Total_N"}, + # "COD": {"method": "_COD"}, + # } + # ) + obj.add_default_units( + { + "time": pyo.units.s, + "length": pyo.units.m, + "mass": pyo.units.kg, + "amount": pyo.units.kmol, + "temperature": pyo.units.K, + } + ) + + +class ASM3PropertiesScaler(CustomScalerBase): + """ + Scaler for the Activated Sludge Model No.1 property package. + + Flow and temperature are scaled by the default value (if no user input provided), and + pressure is scaled assuming an order of magnitude of 1e5 Pa. + """ + + UNIT_SCALING_FACTORS = { + # "QuantityName: (reference units, scaling factor) + "Pressure": (pyo.units.Pa, 1e-6), + } + + DEFAULT_SCALING_FACTORS = { + "flow_vol": 1e1, + "temperature": 1e-1, + } + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + self.scale_variable_by_default(model.temperature, overwrite=overwrite) + self.scale_variable_by_default(model.flow_vol, overwrite=overwrite) + self.scale_variable_by_units(model.pressure, overwrite=overwrite) + + # There are currently no constraints in this model + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + pass + + +class _ASM3StateBlock(StateBlock): + """ + This Class contains methods which should be applied to Property Blocks as a + whole, rather than individual elements of indexed Property Blocks. + """ + + default_scaler = ASM3PropertiesScaler + + def initialize( + self, + state_args=None, + state_vars_fixed=False, + hold_state=False, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + ): + """ + Initialization routine for property package. + + Keyword Arguments: + state_args : Dictionary with initial guesses for the state vars + chosen. Note that if this method is triggered + through the control volume, and if initial guesses + were not provided at the unit model level, the + control volume passes the inlet values as initial + guess.The keys for the state_args dictionary are: + flow_vol : value at which to initialize total volumetric flow (default=None) + alkalinity: value of alkalinity expressed as molar concentration + conc_mass_comp : value at which to initialize component concentrations (default=None) + pressure : value at which to initialize pressure (default=None) + temperature : value at which to initialize temperature (default=None) + outlvl : sets output level of initialization routine + state_vars_fixed: Flag to denote if state vars have already been fixed. + True - states have already been fixed and + initialization does not need to worry + about fixing and unfixing variables. + False - states have not been fixed. The state + block will deal with fixing/unfixing. + optarg : solver options dictionary object (default=None, use + default solver options) + solver : str indicating which solver to use during + initialization (default = None, use default solver) + hold_state : flag indicating whether the initialization routine + should unfix any state variables fixed during + initialization (default=False). + True - states variables are not unfixed, and + a dict of returned containing flags for + which states were fixed during initialization. + False - state variables are unfixed after + initialization by calling the release_state method. + + Returns: + If hold_states is True, returns a dict containing flags for + which states were fixed during initialization. + """ + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="properties") + + if state_vars_fixed is False: + # Fix state variables if not already fixed + flags = fix_state_vars(self, state_args) + + else: + # Check when the state vars are fixed already result in dof 0 + for k in self.keys(): + if degrees_of_freedom(self[k]) != 0: + raise Exception( + "State vars fixed but degrees of freedom " + "for state block is not zero during " + "initialization." + ) + + if state_vars_fixed is False: + if hold_state is True: + return flags + else: + self.release_state(flags) + + init_log.info("Initialization Complete.") + + def release_state(self, flags, outlvl=idaeslog.NOTSET): + """ + Method to relase state variables fixed during initialization. + + Keyword Arguments: + flags : dict containing information of which state variables + were fixed during initialization, and should now be + unfixed. This dict is returned by initialize if + hold_state=True. + outlvl : sets output level of of logging + """ + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="properties") + + if flags is None: + return + # Unfix state variables + revert_state_vars(self, flags) + init_log.info("State Released.") + + +@declare_process_block_class("ASM3StateBlock", block_class=_ASM3StateBlock) +class ASM3StateBlockData(StateBlockData): + """ + StateBlock for calculating thermophysical properties associated with the ASM3 + reaction system. + """ + + def build(self): + """ + Callable method for Block construction + """ + super().build() + + # Create state variables + self.flow_vol = pyo.Var( + initialize=1.0, + domain=pyo.NonNegativeReals, + doc="Total volumetric flowrate", + units=pyo.units.m**3 / pyo.units.s, + ) + self.pressure = pyo.Var( + domain=pyo.NonNegativeReals, + initialize=101325.0, + bounds=(1e3, 1e6), + doc="Pressure", + units=pyo.units.Pa, + ) + self.temperature = pyo.Var( + domain=pyo.NonNegativeReals, + initialize=298.15, + bounds=(273.15, 323.15), + doc="Temperature", + units=pyo.units.K, + ) + self.conc_mass_comp = pyo.Var( + self.params.solute_set, + domain=pyo.NonNegativeReals, + initialize=0.1, + doc="Component mass concentrations", + units=pyo.units.kg / pyo.units.m**3, + ) + self.alkalinity = pyo.Var( + domain=pyo.NonNegativeReals, + initialize=1, + doc="Alkalinity in molar concentration", + units=pyo.units.kmol / pyo.units.m**3, + ) + + # Material and energy flow and density expressions + def material_flow_expression(self, j): + if j == "H2O": + return self.flow_vol * self.params.dens_mass + elif j == "S_ALK": + # Convert moles of alkalinity to mass of C assuming all is HCO3- + return ( + self.flow_vol + * self.alkalinity + * (12 * pyo.units.kg / pyo.units.kmol) + ) + else: + return self.flow_vol * self.conc_mass_comp[j] + + self.material_flow_expression = pyo.Expression( + self.component_list, + rule=material_flow_expression, + doc="Material flow terms", + ) + + def enthalpy_flow_expression(self): + return ( + self.flow_vol + * self.params.dens_mass + * self.params.cp_mass + * (self.temperature - self.params.temperature_ref) + ) + + self.enthalpy_flow_expression = pyo.Expression( + rule=enthalpy_flow_expression, doc="Enthalpy flow term" + ) + + def material_density_expression(self, j): + if j == "H2O": + return self.params.dens_mass + elif j == "S_ALK": + # Convert moles of alkalinity to mass of C assuming all is HCO3- + return self.alkalinity * (12 * pyo.units.kg / pyo.units.kmol) + else: + return self.conc_mass_comp[j] + + self.material_density_expression = pyo.Expression( + self.component_list, + rule=material_density_expression, + doc="Material density terms", + ) + + def energy_density_expression(self): + return ( + self.params.dens_mass + * self.params.cp_mass + * (self.temperature - self.params.temperature_ref) + ) + + self.energy_density_expression = pyo.Expression( + rule=energy_density_expression, doc="Energy density term" + ) + + # def _TSS(self): + # tss = ( + # self.conc_mass_comp["X_S"] + # + self.conc_mass_comp["X_I"] + # + self.conc_mass_comp["X_BH"] + # + self.conc_mass_comp["X_BA"] + # + self.conc_mass_comp["X_P"] + # ) + # return self.params.COD_to_SS * tss + # + # self.TSS = pyo.Expression( + # rule=_TSS, + # doc="Total suspended solids (TSS)", + # ) + # + # def _BOD5(self, i): + # bod5 = ( + # self.conc_mass_comp["S_S"] + # + self.conc_mass_comp["X_S"] + # + (1 - self.params.f_p) + # * (self.conc_mass_comp["X_BH"] + self.conc_mass_comp["X_BA"]) + # ) + # # TODO: 0.25 should be a parameter instead as it changes by influent/effluent + # return self.params.BOD5_factor[i] * bod5 + # + # self.BOD5 = pyo.Expression( + # ["raw", "effluent"], + # rule=_BOD5, + # doc="Five-day Biological Oxygen Demand (BOD5)", + # ) + # + # def _COD(self): + # cod = ( + # self.conc_mass_comp["S_S"] + # + self.conc_mass_comp["S_I"] + # + self.conc_mass_comp["X_S"] + # + self.conc_mass_comp["X_I"] + # + self.conc_mass_comp["X_BH"] + # + self.conc_mass_comp["X_BA"] + # + self.conc_mass_comp["X_P"] + # ) + # return cod + # + # self.COD = pyo.Expression( + # rule=_COD, + # doc="Chemical Oxygen Demand", + # ) + # + # def _TKN(self): + # tkn = ( + # self.conc_mass_comp["S_NH"] + # + self.conc_mass_comp["S_ND"] + # + self.conc_mass_comp["X_ND"] + # + self.params.i_xb + # * (self.conc_mass_comp["X_BH"] + self.conc_mass_comp["X_BA"]) + # + self.params.i_xp + # * (self.conc_mass_comp["X_P"] + self.conc_mass_comp["X_I"]) + # ) + # return tkn + # + # self.TKN = pyo.Expression( + # rule=_TKN, + # doc="Total Kjeldahl Nitrogen", + # ) + # + # def _Total_N(self): + # totaln = self.TKN + self.conc_mass_comp["S_NOX"] + # return totaln + # + # self.Total_N = pyo.Expression( + # rule=_Total_N, + # doc="Total Nitrogen", + # ) + + def get_material_flow_terms(self, p, j): + return self.material_flow_expression[j] + + def get_enthalpy_flow_terms(self, p): + return self.enthalpy_flow_expression + + def get_material_density_terms(self, p, j): + return self.material_density_expression[j] + + def get_energy_density_terms(self, p): + return self.energy_density_expression + + def default_material_balance_type(self): + return MaterialBalanceType.componentPhase + + def default_energy_balance_type(self): + return EnergyBalanceType.enthalpyTotal + + def define_state_vars(self): + return { + "flow_vol": self.flow_vol, + "alkalinity": self.alkalinity, + "conc_mass_comp": self.conc_mass_comp, + "temperature": self.temperature, + "pressure": self.pressure, + } + + def define_display_vars(self): + return { + "Volumetric Flowrate": self.flow_vol, + "Molar Alkalinity": self.alkalinity, + "Mass Concentration": self.conc_mass_comp, + "Temperature": self.temperature, + "Pressure": self.pressure, + } + + def get_material_flow_basis(self): + return MaterialFlowBasis.mass + + def calculate_scaling_factors(self): + # Get default scale factors and do calculations from base classes + super().calculate_scaling_factors() + + # No constraints in this model as yet, just need to set scaling factors + # for expressions + sf_F = iscale.get_scaling_factor(self.flow_vol, default=1e2, warning=True) + sf_T = iscale.get_scaling_factor(self.temperature, default=1e-2, warning=True) + + # Mass flow and density terms + for j in self.component_list: + if j == "H2O": + sf_C = pyo.value(1 / self.params.dens_mass) + elif j == "S_ALK": + sf_C = 1e-1 * iscale.get_scaling_factor( + self.alkalinity, default=1, warning=True + ) + else: + sf_C = iscale.get_scaling_factor( + self.conc_mass_comp[j], default=1e2, warning=True + ) + + iscale.set_scaling_factor(self.material_flow_expression[j], sf_F * sf_C) + iscale.set_scaling_factor(self.material_density_expression[j], sf_C) + + # Enthalpy and energy terms + sf_rho_cp = pyo.value(1 / (self.params.dens_mass * self.params.cp_mass)) + iscale.set_scaling_factor( + self.enthalpy_flow_expression, sf_F * sf_rho_cp * sf_T + ) + iscale.set_scaling_factor(self.energy_density_expression, sf_rho_cp * sf_T) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py new file mode 100644 index 0000000000..7c0fd99d1f --- /dev/null +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py @@ -0,0 +1,603 @@ +################################################################################# +# 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/" +################################################################################# +""" +ASM1 reaction package. + +References: + +[1] Henze, M., Grady, C.P.L., Gujer, W., Marais, G.v.R., Matsuo, T., +"Activated Sludge Model No. 1", 1987, IAWPRC Task Group on Mathematical Modeling +for Design and Operation of Biological Wastewater Treatment +[2] J. Alex, L. Benedetti, J. Copp, K.V. Gernaey, U. Jeppsson, I. Nopens, M.N. Pons, +J.P. Steyer and P. Vanrolleghem, "Benchmark Simulation Model no. 1 (BSM1)", 2018 +""" + +# Import Pyomo libraries +import pyomo.environ as pyo + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + MaterialFlowBasis, + ReactionParameterBlock, + ReactionBlockDataBase, + ReactionBlockBase, +) +from idaes.core.util.misc import add_object_reference +from idaes.core.util.exceptions import BurntToast +import idaes.logger as idaeslog +from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme +import idaes.core.util.scaling as iscale + +# Some more information about this module +__author__ = "Chenyu Wang, Adam Atia" + + +# Set up logger +_log = idaeslog.getLogger(__name__) + + +@declare_process_block_class("ASM3ReactionParameterBlock") +class ASM3ReactionParameterData(ReactionParameterBlock): + """ + Reaction Parameter Block Class + """ + + def build(self): + """ + Callable method for Block construction. + """ + super().build() + + self._reaction_block_class = ASM3ReactionBlock + + # Reaction Index + # Reaction names based on standard numbering in ASM3 paper + + # R1: Hydrolysis + + # Heterotrophic organisms, aerobic and denitrifying activity + # R2: Aerobic storage of S_S + # R3: Anoxic storage of S_S + # R4: Aerobic growth of X_H + # R5: Anoxic growth (denitrific) + # R6: Aerobic endog. respiration + # R7: Anoxic endog. respiration + # R8: Aerobic respiration of X_STO + # R9: Anoxic respiration of X_STO + + # Autotrophic organisms, nitrifying activity + # R10: Aerobic growth of X_A + # R11: Aerobic endog. respiration + # R12: Anoxic endog. respiration + + self.rate_reaction_idx = pyo.Set( + initialize=[ + "R1", + "R2", + "R3", + "R4", + "R5", + "R6", + "R7", + "R8", + "R9", + "R10", + "R10", + "R11", + "R12", + ] + ) + + # Stoichiometric Parameters + self.Y_A = pyo.Var( + initialize=0.24, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Yield of autotrophic biomass per N03-N (g-COD-X_A / g-N-S_NOX)", + ) + self.Y_H_O2 = pyo.Var( + initialize=0.63, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Aerobic yield of heterotrophic biomass (g-COD-X_H / g-COD-X_STO)", + ) + self.Y_H_NOX = pyo.Var( + initialize=0.54, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Anoxic yield of heterotrophic biomass (g-COD-X_H / g-COD-X_STO)", + ) + self.Y_STO_O2 = pyo.Var( + initialize=0.85, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Aerobic yield of stored product per S_S (g-COD-X_STO / g-COD-S_S)", + ) + self.Y_STO_NOX = pyo.Var( + initialize=0.80, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Anoxic yield of stored product per per S_S (g-COD-X_STO / g-COD-S_S)", + ) + + add_object_reference(self, "f_SI", self.config.property_package.f_SI) + add_object_reference(self, "i_NSI", self.config.property_package.i_NSI) + add_object_reference(self, "i_NSS", self.config.property_package.i_NSS) + add_object_reference(self, "i_NXI", self.config.property_package.i_NXI) + add_object_reference(self, "i_NXS", self.config.property_package.i_NXS) + add_object_reference(self, "i_NBM", self.config.property_package.i_NBM) + add_object_reference(self, "i_TSXI", self.config.property_package.i_TSXI) + add_object_reference(self, "i_TSXS", self.config.property_package.i_TSXS) + add_object_reference(self, "i_TSBM", self.config.property_package.i_TSBM) + add_object_reference(self, "i_TSSTO", self.config.property_package.i_TSSTO) + + # Kinetic Parameters + k_H_dict = {"10C": 2, "20C": 3} + self.k_H = pyo.Var( + k_H_dict.keys(), + domain=pyo.PositiveReals, + initialize=k_H_dict, + units=pyo.units.day**-1, + doc="Hydrolysis rate constant (g-COD-X_S / g-COD-X_H / day)", + ) + + self.mu_A = pyo.Var( + initialize=0.5, + units=1 / pyo.units.day, + domain=pyo.PositiveReals, + doc="Maximum specific growth rate for autotrophic biomass, mu_A", + ) + self.mu_H = pyo.Var( + initialize=4, + units=1 / pyo.units.day, + domain=pyo.PositiveReals, + doc="Maximum specific growth rate for heterotrophic biomass, mu_H", + ) + self.K_S = pyo.Var( + initialize=10e-3, + units=pyo.units.kg / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Half-saturation coefficient for heterotrophic biomass, K_S", + ) + self.K_OH = pyo.Var( + initialize=0.2e-3, + units=pyo.units.kg / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Oxygen half-saturation coefficient for heterotrophic biomass, K_O,H", + ) + self.K_OA = pyo.Var( + initialize=0.4e-3, + units=pyo.units.kg / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Oxygen half-saturation coefficient for autotrophic biomass, K_O,A", + ) + self.K_NO = pyo.Var( + initialize=0.5e-3, + units=pyo.units.kg / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Nitrate half-saturation coefficient for denitrifying heterotrophic biomass, K_NO", + ) + self.b_H = pyo.Var( + initialize=0.3, + units=1 / pyo.units.day, + domain=pyo.PositiveReals, + doc="Decay coefficient for heterotrophic biomass, b_H", + ) + self.b_A = pyo.Var( + initialize=0.05, + units=1 / pyo.units.day, + domain=pyo.PositiveReals, + doc="Decay coefficient for autotrophic biomass, b_A", + ) + self.eta_g = pyo.Var( + initialize=0.8, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Correction factor for mu_H under anoxic conditions, eta_g", + ) + self.eta_h = pyo.Var( + initialize=0.8, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Correction factor for hydrolysis under anoxic conditions, eta_h", + ) + self.k_h = pyo.Var( + initialize=3, + units=1 / pyo.units.day, + domain=pyo.PositiveReals, + doc="Maximum specific hydrolysis rate, k_h", + ) + self.K_X = pyo.Var( + initialize=0.1, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Half-saturation coefficient for hydrolysis of slowly biodegradable substrate, K_X", + ) + self.K_NH = pyo.Var( + initialize=1e-3, + units=pyo.units.kg / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Ammonia Half-saturation coefficient for autotrophic biomass, K_NH", + ) + self.k_a = pyo.Var( + initialize=0.05e3, + units=pyo.units.m**3 / pyo.units.kg / pyo.units.day, + domain=pyo.PositiveReals, + doc="Ammonification rate, k_a", + ) + + # Reaction Stoichiometry + # This is the stoichiometric part the Peterson matrix in dict form + # Note that reaction stoichiometry is on a mass basis. + # For alkalinity, this requires converting the mass of nitrogen species + # reacted to mass of alkalinity converted using a charge balance (effectively MW_C/MW_N) + mw_alk = 12 * pyo.units.kg / pyo.units.kmol + mw_n = 14 * pyo.units.kg / pyo.units.kmol + self.rate_reaction_stoichiometry = { + # R1: Aerobic growth of heterotrophs + ("R1", "Liq", "H2O"): 0, + ("R1", "Liq", "S_I"): 0, + ("R1", "Liq", "S_S"): -1 / self.Y_H, + ("R1", "Liq", "X_I"): 0, + ("R1", "Liq", "X_S"): 0, + ("R1", "Liq", "X_BH"): 1, + ("R1", "Liq", "X_BA"): 0, + ("R1", "Liq", "X_P"): 0, + ("R1", "Liq", "S_O"): -(1 - self.Y_H) / self.Y_H, + ("R1", "Liq", "S_NO"): 0, + ("R1", "Liq", "S_NH"): -self.i_xb, + ("R1", "Liq", "S_ND"): 0, + ("R1", "Liq", "X_ND"): 0, + ("R1", "Liq", "S_ALK"): -self.i_xb * (mw_alk / mw_n), + # R2: Anoxic growth of heterotrophs + ("R2", "Liq", "H2O"): 0, + ("R2", "Liq", "S_I"): 0, + ("R2", "Liq", "S_S"): -1 / self.Y_H, + ("R2", "Liq", "X_I"): 0, + ("R2", "Liq", "X_S"): 0, + ("R2", "Liq", "X_BH"): 1, + ("R2", "Liq", "X_BA"): 0, + ("R2", "Liq", "X_P"): 0, + ("R2", "Liq", "S_O"): 0, + ("R2", "Liq", "S_NO"): -(1 - self.Y_H) / (2.86 * self.Y_H), + ("R2", "Liq", "S_NH"): -self.i_xb, + ("R2", "Liq", "S_ND"): 0, + ("R2", "Liq", "X_ND"): 0, + ("R2", "Liq", "S_ALK"): ((1 - self.Y_H) / (2.86 * self.Y_H) - self.i_xb) + * (mw_alk / mw_n), + # R3: Aerobic growth of autotrophs + ("R3", "Liq", "H2O"): 0, + ("R3", "Liq", "S_I"): 0, + ("R3", "Liq", "S_S"): 0, + ("R3", "Liq", "X_I"): 0, + ("R3", "Liq", "X_S"): 0, + ("R3", "Liq", "X_BH"): 0, + ("R3", "Liq", "X_BA"): 1, + ("R3", "Liq", "X_P"): 0, + ("R3", "Liq", "S_O"): -(4.57 - self.Y_A) / self.Y_A, + ("R3", "Liq", "S_NO"): 1 / self.Y_A, + ("R3", "Liq", "S_NH"): -self.i_xb - 1 / self.Y_A, + ("R3", "Liq", "S_ND"): 0, + ("R3", "Liq", "X_ND"): 0, + ("R3", "Liq", "S_ALK"): (-self.i_xb - 2 / (self.Y_A)) * (mw_alk / mw_n), + # R4: Decay of heterotrophs + ("R4", "Liq", "H2O"): 0, + ("R4", "Liq", "S_I"): 0, + ("R4", "Liq", "S_S"): 0, + ("R4", "Liq", "X_I"): 0, + ("R4", "Liq", "X_S"): 1 - self.f_p, + ("R4", "Liq", "X_BH"): -1, + ("R4", "Liq", "X_BA"): 0, + ("R4", "Liq", "X_P"): self.f_p, + ("R4", "Liq", "S_O"): 0, + ("R4", "Liq", "S_NO"): 0, + ("R4", "Liq", "S_NH"): 0, + ("R4", "Liq", "S_ND"): 0, + ("R4", "Liq", "X_ND"): self.i_xb - self.f_p * self.i_xp, + ("R4", "Liq", "S_ALK"): 0, + # R5: Decay of autotrophs + ("R5", "Liq", "H2O"): 0, + ("R5", "Liq", "S_I"): 0, + ("R5", "Liq", "S_S"): 0, + ("R5", "Liq", "X_I"): 0, + ("R5", "Liq", "X_S"): 1 - self.f_p, + ("R5", "Liq", "X_BH"): 0, + ("R5", "Liq", "X_BA"): -1, + ("R5", "Liq", "X_P"): self.f_p, + ("R5", "Liq", "S_O"): 0, + ("R5", "Liq", "S_NO"): 0, + ("R5", "Liq", "S_NH"): 0, + ("R5", "Liq", "S_ND"): 0, + ("R5", "Liq", "X_ND"): self.i_xb - self.f_p * self.i_xp, + ("R5", "Liq", "S_ALK"): 0, + # R6: Ammonification of soluble organic nitrogen + ("R6", "Liq", "H2O"): 0, + ("R6", "Liq", "S_I"): 0, + ("R6", "Liq", "S_S"): 0, + ("R6", "Liq", "X_I"): 0, + ("R6", "Liq", "X_S"): 0, + ("R6", "Liq", "X_BH"): 0, + ("R6", "Liq", "X_BA"): 0, + ("R6", "Liq", "X_P"): 0, + ("R6", "Liq", "S_O"): 0, + ("R6", "Liq", "S_NO"): 0, + ("R6", "Liq", "S_NH"): 1, + ("R6", "Liq", "S_ND"): -1, + ("R6", "Liq", "X_ND"): 0, + ("R6", "Liq", "S_ALK"): 1 * (mw_alk / mw_n), + # R7: Hydrolysis of entrapped organics + ("R7", "Liq", "H2O"): 0, + ("R7", "Liq", "S_I"): 0, + ("R7", "Liq", "S_S"): 1, + ("R7", "Liq", "X_I"): 0, + ("R7", "Liq", "X_S"): -1, + ("R7", "Liq", "X_BH"): 0, + ("R7", "Liq", "X_BA"): 0, + ("R7", "Liq", "X_P"): 0, + ("R7", "Liq", "S_O"): 0, + ("R7", "Liq", "S_NO"): 0, + ("R7", "Liq", "S_NH"): 0, + ("R7", "Liq", "S_ND"): 0, + ("R7", "Liq", "X_ND"): 0, + ("R7", "Liq", "S_ALK"): 0, + # R8: Hydrolysis of entrapped organic nitrogen + ("R8", "Liq", "H2O"): 0, + ("R8", "Liq", "S_I"): 0, + ("R8", "Liq", "S_S"): 0, + ("R8", "Liq", "X_I"): 0, + ("R8", "Liq", "X_S"): 0, + ("R8", "Liq", "X_BH"): 0, + ("R8", "Liq", "X_BA"): 0, + ("R8", "Liq", "X_P"): 0, + ("R8", "Liq", "S_O"): 0, + ("R8", "Liq", "S_NO"): 0, + ("R8", "Liq", "S_NH"): 0, + ("R8", "Liq", "S_ND"): 1, + ("R8", "Liq", "X_ND"): -1, + ("R8", "Liq", "S_ALK"): 0, + } + # Fix all the variables we just created + for v in self.component_objects(pyo.Var, descend_into=False): + v.fix() + + @classmethod + def define_metadata(cls, obj): + obj.add_properties( + { + "reaction_rate": {"method": "_rxn_rate"}, + } + ) + obj.add_default_units( + { + "time": pyo.units.s, + "length": pyo.units.m, + "mass": pyo.units.kg, + "amount": pyo.units.kmol, + "temperature": pyo.units.K, + } + ) + + +class ASM3ReactionScaler(CustomScalerBase): + """ + Scaler for the Activated Sludge Model No.1 reaction package. + + Variables are scaled by their default scaling factor (if no user input provided), and constraints + are scaled using the inverse maximum scheme. + """ + + # TODO: Revisit this scaling factor + DEFAULT_SCALING_FACTORS = {"reaction_rate": 1e2} + + def variable_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + + if model.is_property_constructed("reaction_rate"): + for j in model.reaction_rate.values(): + self.scale_variable_by_default(j, overwrite=overwrite) + + def constraint_scaling_routine( + self, model, overwrite: bool = False, submodel_scalers: dict = None + ): + # TODO: Revisit this scaling methodology + # Consider scale_constraint_by_default or scale_constraints_by_jacobian_norm + if model.is_property_constructed("rate_expression"): + for j in model.rate_expression.values(): + self.scale_constraint_by_nominal_value( + j, + scheme=ConstraintScalingScheme.inverseMaximum, + overwrite=overwrite, + ) + + +class _ASM3ReactionBlock(ReactionBlockBase): + """ + This Class contains methods which should be applied to Reaction Blocks as a + whole, rather than individual elements of indexed Reaction Blocks. + """ + + default_scaler = ASM3ReactionScaler + + def initialize(self, outlvl=idaeslog.NOTSET, **kwargs): + """ + Initialization routine for reaction package. + + Keyword Arguments: + outlvl : sets output level of initialization routine + + Returns: + None + """ + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="properties") + init_log.info("Initialization Complete.") + + +@declare_process_block_class("ASM3ReactionBlock", block_class=_ASM3ReactionBlock) +class ASM3ReactionBlockData(ReactionBlockDataBase): + """ + ReactionBlock for ASM3. + """ + + def build(self): + """ + Callable method for Block construction + """ + super().build() + + # Create references to state vars + # Concentration + add_object_reference(self, "conc_mass_comp_ref", self.state_ref.conc_mass_comp) + + # Rate of reaction method + def _rxn_rate(self): + self.reaction_rate = pyo.Var( + self.params.rate_reaction_idx, + initialize=0, + doc="Rate of reaction", + units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + + try: + + def rate_expression_rule(b, r): + if r == "R1": + # R1: Aerobic growth of heterotrophs + return b.reaction_rate[r] == pyo.units.convert( + b.params.mu_H + * ( + b.conc_mass_comp_ref["S_S"] + / (b.params.K_S + b.conc_mass_comp_ref["S_S"]) + ) + * ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) + ) + * b.conc_mass_comp_ref["X_BH"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R2": + # R2: Anoxic growth of heterotrophs + return b.reaction_rate[r] == pyo.units.convert( + b.params.mu_H + * ( + b.conc_mass_comp_ref["S_S"] + / (b.params.K_S + b.conc_mass_comp_ref["S_S"]) + ) + * ( + b.params.K_OH + / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) + ) + * ( + b.conc_mass_comp_ref["S_NO"] + / (b.params.K_NO + b.conc_mass_comp_ref["S_NO"]) + ) + * b.params.eta_g + * b.conc_mass_comp_ref["X_BH"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R3": + # R3: Aerobic growth of autotrophs + return b.reaction_rate[r] == pyo.units.convert( + b.params.mu_A + * ( + b.conc_mass_comp_ref["S_NH"] + / (b.params.K_NH + b.conc_mass_comp_ref["S_NH"]) + ) + * ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_OA + b.conc_mass_comp_ref["S_O"]) + ) + * b.conc_mass_comp_ref["X_BA"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R4": + # R4: Decay of heterotrophs + return b.reaction_rate[r] == pyo.units.convert( + b.params.b_H * b.conc_mass_comp_ref["X_BH"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R5": + # R5: Decay of autotrophs + return b.reaction_rate[r] == pyo.units.convert( + b.params.b_A * b.conc_mass_comp_ref["X_BA"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R6": + # R6: Ammonification of soluble organic nitrogen + return b.reaction_rate[r] == pyo.units.convert( + b.params.k_a + * b.conc_mass_comp_ref["S_ND"] + * b.conc_mass_comp_ref["X_BH"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R7": + # 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.params.K_X * b.conc_mass_comp_ref["X_BH"] + + b.conc_mass_comp_ref["X_S"] + ) + * ( + ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) + ) + + b.params.eta_h + * b.params.K_OH + / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) + * ( + b.conc_mass_comp_ref["S_NO"] + / (b.params.K_NO + b.conc_mass_comp_ref["S_NO"]) + ) + ) + * b.conc_mass_comp_ref["X_BH"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R8": + # 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"] + + 1e-10 * pyo.units.kg / pyo.units.m**3 + ) + ) + ) + else: + raise BurntToast() + + self.rate_expression = pyo.Constraint( + self.params.rate_reaction_idx, + rule=rate_expression_rule, + doc="ASM3 rate expressions", + ) + + except AttributeError: + # If constraint fails, clean up so that DAE can try again later + self.del_component(self.reaction_rate) + self.del_component(self.rate_expression) + raise + + def get_reaction_rate_basis(self): + return MaterialFlowBasis.mass + + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + iscale.constraint_scaling_transform(self.rate_expression["R5"], 1e3) + iscale.constraint_scaling_transform(self.rate_expression["R3"], 1e3) + iscale.constraint_scaling_transform(self.rate_expression["R4"], 1e3) From cb4c637c51e8b445eda32b8ce593cf24f97ae022 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Wed, 11 Jun 2025 17:26:38 -0400 Subject: [PATCH 07/21] Completed reaction package and properties package --- .../activated_sludge/asm3_properties.py | 22 +- .../activated_sludge/asm3_reactions.py | 691 ++++++++++++------ 2 files changed, 490 insertions(+), 223 deletions(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py index 5b403f4ea2..79521215ac 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py @@ -67,8 +67,8 @@ def build(self): self.S_O = Solute(doc="Dissolved Oxygen, S_O") self.S_I = Solute(doc="Inert soluble organic material, S_I") self.S_S = Solute(doc="Readily biodegradable substrates, S_S") - self.S_NH = Solute( - doc="Ammonium plus ammonia nitrogen (NH4+ + NH3 nitrogen), S_NH" + self.S_NH4 = Solute( + doc="Ammonium plus ammonia nitrogen (NH4+ + NH3 nitrogen), S_NH4" ) self.S_N2 = Solute(doc="Dinitrogen (N2), S_N2") self.S_NOX = Solute( @@ -92,7 +92,7 @@ def build(self): "S_O", "S_I", "S_S", - "S_NH", + "S_NH4", "S_N2", "S_NOX", "S_ALK", @@ -139,6 +139,12 @@ def build(self): domain=pyo.NonNegativeReals, doc="Production of S_I in hydrolysis (g-COD-S_I / g-COD-X_S)", ) + self.f_XI = pyo.Var( + initialize=0.2, + units=pyo.units.dimensionless, + domain=pyo.NonNegativeReals, + doc="Production of X_I in endog. respiration (g-COD-X_I / g-COD-X_BM)", + ) self.i_NSI = pyo.Var( initialize=0.01, units=pyo.units.dimensionless, @@ -169,25 +175,25 @@ def build(self): domain=pyo.PositiveReals, doc="N content of biomass, X_H, X_A (g-N / g-COD-X_H or X_A)", ) - self.i_TSXI = pyo.Var( + self.i_SSXI = pyo.Var( initialize=0.75, units=pyo.units.dimensionless, domain=pyo.PositiveReals, doc="TSS to COD ratio for X_I (g-TSS / g-COD-X_I)", ) - self.i_TSXS = pyo.Var( + self.i_SSXS = pyo.Var( initialize=0.75, units=pyo.units.dimensionless, domain=pyo.PositiveReals, doc="TSS to COD ratio for X_S (g-TSS / g-COD-X_S)", ) - self.i_TSBM = pyo.Var( + self.i_SSBM = pyo.Var( initialize=0.90, units=pyo.units.dimensionless, domain=pyo.PositiveReals, doc="TSS to COD ratio for for biomass,X_H, X_A (g-TSS / g-COD-X_H or X_A)", ) - self.i_TSSTO = pyo.Var( + self.i_SSSTO = pyo.Var( initialize=0.60, units=pyo.units.dimensionless, domain=pyo.PositiveReals, @@ -547,7 +553,7 @@ def energy_density_expression(self): # # def _TKN(self): # tkn = ( - # self.conc_mass_comp["S_NH"] + # self.conc_mass_comp["S_NH4"] # + self.conc_mass_comp["S_ND"] # + self.conc_mass_comp["X_ND"] # + self.params.i_xb diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py index 7c0fd99d1f..2c536c5af8 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py @@ -131,15 +131,16 @@ def build(self): ) add_object_reference(self, "f_SI", self.config.property_package.f_SI) + add_object_reference(self, "f_XI", self.config.property_package.f_XI) add_object_reference(self, "i_NSI", self.config.property_package.i_NSI) add_object_reference(self, "i_NSS", self.config.property_package.i_NSS) add_object_reference(self, "i_NXI", self.config.property_package.i_NXI) add_object_reference(self, "i_NXS", self.config.property_package.i_NXS) add_object_reference(self, "i_NBM", self.config.property_package.i_NBM) - add_object_reference(self, "i_TSXI", self.config.property_package.i_TSXI) - add_object_reference(self, "i_TSXS", self.config.property_package.i_TSXS) - add_object_reference(self, "i_TSBM", self.config.property_package.i_TSBM) - add_object_reference(self, "i_TSSTO", self.config.property_package.i_TSSTO) + add_object_reference(self, "i_SSXI", self.config.property_package.i_SSXI) + add_object_reference(self, "i_SSXS", self.config.property_package.i_SSXS) + add_object_reference(self, "i_SSBM", self.config.property_package.i_SSBM) + add_object_reference(self, "i_SSSTO", self.config.property_package.i_SSSTO) # Kinetic Parameters k_H_dict = {"10C": 2, "20C": 3} @@ -150,92 +151,200 @@ def build(self): units=pyo.units.day**-1, doc="Hydrolysis rate constant (g-COD-X_S / g-COD-X_H / day)", ) + self.K_X = pyo.Var( + initialize=1, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Hydrolysis saturation constant (g-COD-X_S / g-COD-X_H)", + ) - self.mu_A = pyo.Var( + # Heterotrophic organisms X_H, aerobic and denitrifying activity + k_STO_dict = {"10C": 2.5, "20C": 5} + self.k_STO = pyo.Var( + k_STO_dict.keys(), + domain=pyo.PositiveReals, + initialize=k_STO_dict, + units=pyo.units.day**-1, + doc="Storage rate constant (g-COD-S_S / g-COD-X_H / day)", + ) + self.eta_NOX = pyo.Var( + initialize=0.6, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Anoxic reduction factor", + ) + self.K_O2 = pyo.Var( + initialize=0.2, + units=pyo.units.g / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Saturation constant for S_NO2 (g-O2 / m3)", + ) + self.K_NOX = pyo.Var( initialize=0.5, - units=1 / pyo.units.day, + units=pyo.units.g / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Saturation constant for S_NOX (g-N-NO3 / m3)", + ) + self.K_S = pyo.Var( + initialize=2, + units=pyo.units.g / pyo.units.m**3, + domain=pyo.PositiveReals, + doc="Saturation constant for for substrate S_S (g-COD-S_S / m3)", + ) + self.K_STO = pyo.Var( + initialize=1, + units=pyo.units.dimensionless, domain=pyo.PositiveReals, - doc="Maximum specific growth rate for autotrophic biomass, mu_A", + doc="Saturation constant for for X_STO (g-COD-X_STO / g-COD-X_H)", ) + mu_H_dict = {"10C": 1, "20C": 2} self.mu_H = pyo.Var( - initialize=4, - units=1 / pyo.units.day, + mu_H_dict.keys(), domain=pyo.PositiveReals, - doc="Maximum specific growth rate for heterotrophic biomass, mu_H", + initialize=mu_H_dict, + units=pyo.units.day**-1, + doc="Heterotrophic max. growth rate of X_H (day^-1)", ) - self.K_S = pyo.Var( - initialize=10e-3, - units=pyo.units.kg / pyo.units.m**3, + self.K_NH4 = pyo.Var( + initialize=0.01, + units=pyo.units.g / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Half-saturation coefficient for heterotrophic biomass, K_S", + doc="Saturation constant for ammonium, S_NH4 (g-N / m3)", ) - self.K_OH = pyo.Var( - initialize=0.2e-3, - units=pyo.units.kg / pyo.units.m**3, + self.K_ALK = pyo.Var( + initialize=0.1, + units=pyo.units.mol / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Oxygen half-saturation coefficient for heterotrophic biomass, K_O,H", + doc="Saturation constant for alkalinity for X_H (mol-HCO3- / m3)", ) - self.K_OA = pyo.Var( - initialize=0.4e-3, - units=pyo.units.kg / pyo.units.m**3, + b_H_O2_dict = {"10C": 0.1, "20C": 0.2} + self.b_H_O2 = pyo.Var( + b_H_O2_dict.keys(), domain=pyo.PositiveReals, - doc="Oxygen half-saturation coefficient for autotrophic biomass, K_O,A", + initialize=b_H_O2_dict, + units=pyo.units.day**-1, + doc="Aerobic endogenous respiration rate of X_H (day^-1)", ) - self.K_NO = pyo.Var( - initialize=0.5e-3, - units=pyo.units.kg / pyo.units.m**3, + b_H_NOX_dict = {"10C": 0.05, "20C": 0.1} + self.b_H_NOX = pyo.Var( + b_H_NOX_dict.keys(), domain=pyo.PositiveReals, - doc="Nitrate half-saturation coefficient for denitrifying heterotrophic biomass, K_NO", + initialize=b_H_NOX_dict, + units=pyo.units.day**-1, + doc="Anoxic endogenous respiration rate of X_H (day^-1)", ) - self.b_H = pyo.Var( - initialize=0.3, - units=1 / pyo.units.day, + b_STO_O2_dict = {"10C": 0.1, "20C": 0.2} + self.b_STO_O2 = pyo.Var( + b_STO_O2_dict.keys(), domain=pyo.PositiveReals, - doc="Decay coefficient for heterotrophic biomass, b_H", + initialize=b_STO_O2_dict, + units=pyo.units.day**-1, + doc="Aerobic respiration rate for X_STO (day^-1)", ) - self.b_A = pyo.Var( - initialize=0.05, - units=1 / pyo.units.day, + b_STO_NOX_dict = {"10C": 0.05, "20C": 0.1} + self.b_STO_NOX = pyo.Var( + b_STO_NOX_dict.keys(), domain=pyo.PositiveReals, - doc="Decay coefficient for autotrophic biomass, b_A", + initialize=b_STO_NOX_dict, + units=pyo.units.day**-1, + doc="Anoxic respiration rate for X_STO (day^-1)", ) - self.eta_g = pyo.Var( - initialize=0.8, - units=pyo.units.dimensionless, + + # Autotrophic organisms X_A, nitrifying activity + mu_A_dict = {"10C": 0.35, "20C": 1} + self.mu_A = pyo.Var( + mu_A_dict.keys(), domain=pyo.PositiveReals, - doc="Correction factor for mu_H under anoxic conditions, eta_g", + initialize=mu_A_dict, + units=pyo.units.day**-1, + doc="Autotrophic max. growth rate of X_A (day^-1)", ) - self.eta_h = pyo.Var( - initialize=0.8, - units=pyo.units.dimensionless, + self.K_A_NH4 = pyo.Var( + initialize=1, + units=pyo.units.g / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Correction factor for hydrolysis under anoxic conditions, eta_h", + doc="Ammonium substrate saturation for X_A (g-N / m3)", ) - self.k_h = pyo.Var( - initialize=3, - units=1 / pyo.units.day, + self.K_A_O2 = pyo.Var( + initialize=0.5, + units=pyo.units.g / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Maximum specific hydrolysis rate, k_h", + doc="Oxygen saturation for nitrifiers (g-O2 / m3)", ) - self.K_X = pyo.Var( - initialize=0.1, - units=pyo.units.dimensionless, + self.K_A_ALK = pyo.Var( + initialize=0.5, + units=pyo.units.mol / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Half-saturation coefficient for hydrolysis of slowly biodegradable substrate, K_X", + doc="Bicarbonate saturation for nitrifiers (mol-HCO3- / m3)", ) - self.K_NH = pyo.Var( - initialize=1e-3, - units=pyo.units.kg / pyo.units.m**3, + b_A_O2_dict = {"10C": 0.05, "20C": 0.15} + self.b_A_O2 = pyo.Var( + b_A_O2_dict.keys(), domain=pyo.PositiveReals, - doc="Ammonia Half-saturation coefficient for autotrophic biomass, K_NH", + initialize=b_A_O2_dict, + units=pyo.units.day**-1, + doc="Aerobic endogenous respiration rate of X_A (day^-1)", ) - self.k_a = pyo.Var( - initialize=0.05e3, - units=pyo.units.m**3 / pyo.units.kg / pyo.units.day, + b_A_NOX_dict = {"10C": 0.02, "20C": 0.05} + self.b_A_NOX = pyo.Var( + b_A_NOX_dict.keys(), domain=pyo.PositiveReals, - doc="Ammonification rate, k_a", + initialize=b_A_NOX_dict, + units=pyo.units.day**-1, + doc="Anoxic endogenous respiration rate of X_A (day^-1)", ) + # Stoichiometric numbers from Table 1 + # obtained by \sum_i^12 νji*ikI + x1 = 1.0 - self.f_SI + x2 = -1.0 + self.Y_STO_O2 + x3 = (-1.0 + self.Y_STO_NOX) / (64.0 / 14.0 - 24.0 / 14.0) + x4 = 1.0 - 1.0 / self.Y_H_O2 + x5 = (+1.0 - 1.0 / self.Y_H_NOX) / (64.0 / 14.0 - 24.0 / 14.0) + x6 = -1.0 + self.f_XI + x7 = (self.f_XI - 1.0) / (64.0 / 14.0 - 24.0 / 14.0) + x8 = -1.0 + x9 = -1.0 / (64.0 / 14.0 - 24.0 / 14.0) + x10 = -(64.0 / 14.0) / self.Y_A + 1.0 + x11 = self.f_XI - 1.0 + x12 = (self.f_XI - 1.0) / (64.0 / 14.0 - 24.0 / 14.0) + + y1 = -self.f_SI * self.i_NSI - (1.0 - self.f_SI) * self.i_NSS + self.i_NXS + y2 = self.i_NSS + y3 = self.i_NSS + y4 = -self.i_NBM + y5 = -self.i_NBM + y6 = -self.f_XI * self.i_NXI + self.i_NBM + y7 = -self.f_XI * self.i_NXI + self.i_NBM + y10 = -1.0 / self.Y_A - self.i_NBM + y11 = -self.f_XI * self.i_NXI + self.i_NBM + y12 = -self.f_XI * self.i_NXI + self.i_NBM + + z1 = y1 / 14.0 + z2 = y2 / 14.0 + z3 = y3 / 14.0 - x3 / 14.0 + z4 = y4 / 14.0 + z5 = y5 / 14.0 - x5 / 14.0 + z6 = y6 / 14.0 + z7 = y7 / 14.0 - x7 / 14.0 + z9 = -x9 / 14.0 + z10 = y10 / 14.0 - 1.0 / (self.Y_A * 14.0) + z11 = y11 / 14.0 + z12 = y12 / 14.0 - x12 / 14.0 + + t1 = -self.i_SSXS + t2 = self.Y_STO_O2 * self.i_SSSTO + t3 = self.Y_STO_NOX * self.i_SSSTO + t4 = self.i_SSBM - 1.0 / self.Y_H_O2 * self.i_SSSTO + t5 = self.i_SSBM - 1.0 / self.Y_H_NOX * self.i_SSSTO + t6 = self.f_XI * self.i_SSXI - self.i_SSBM + t7 = self.f_XI * self.i_SSXI - self.i_SSBM + t8 = -self.i_SSSTO + t9 = -self.i_SSSTO + t10 = self.i_SSBM + t11 = self.f_XI * self.i_SSXI - self.i_SSBM + t12 = self.f_XI * self.i_SSXI - self.i_SSBM + # Reaction Stoichiometry # This is the stoichiometric part the Peterson matrix in dict form # Note that reaction stoichiometry is on a mass basis. @@ -244,128 +353,190 @@ def build(self): mw_alk = 12 * pyo.units.kg / pyo.units.kmol mw_n = 14 * pyo.units.kg / pyo.units.kmol self.rate_reaction_stoichiometry = { - # R1: Aerobic growth of heterotrophs + # R1: Hydrolysis ("R1", "Liq", "H2O"): 0, - ("R1", "Liq", "S_I"): 0, - ("R1", "Liq", "S_S"): -1 / self.Y_H, + ("R1", "Liq", "S_O"): 0, + ("R1", "Liq", "S_I"): self.f_SI, + ("R1", "Liq", "S_S"): x1, + ("R1", "Liq", "S_NH4"): y1, + ("R1", "Liq", "S_N2"): 0, + ("R1", "Liq", "S_NOX"): 0, + ("R1", "Liq", "S_ALK"): z1, ("R1", "Liq", "X_I"): 0, - ("R1", "Liq", "X_S"): 0, - ("R1", "Liq", "X_BH"): 1, - ("R1", "Liq", "X_BA"): 0, - ("R1", "Liq", "X_P"): 0, - ("R1", "Liq", "S_O"): -(1 - self.Y_H) / self.Y_H, - ("R1", "Liq", "S_NO"): 0, - ("R1", "Liq", "S_NH"): -self.i_xb, - ("R1", "Liq", "S_ND"): 0, - ("R1", "Liq", "X_ND"): 0, - ("R1", "Liq", "S_ALK"): -self.i_xb * (mw_alk / mw_n), - # R2: Anoxic growth of heterotrophs + ("R1", "Liq", "X_S"): -1, + ("R1", "Liq", "X_H"): 0, + ("R1", "Liq", "X_STO"): 0, + ("R1", "Liq", "X_A"): 0, + ("R1", "Liq", "X_TSS"): -self.i_SSXS, + # Heterotrophic organisms, aerobic and denitrifying activity + # R2: Aerobic storage of S_S ("R2", "Liq", "H2O"): 0, + ("R2", "Liq", "S_O"): x2, ("R2", "Liq", "S_I"): 0, - ("R2", "Liq", "S_S"): -1 / self.Y_H, + ("R2", "Liq", "S_S"): -1, + ("R2", "Liq", "S_NH4"): y2, + ("R2", "Liq", "S_N2"): 0, + ("R2", "Liq", "S_NOX"): 0, + ("R2", "Liq", "S_ALK"): z2, ("R2", "Liq", "X_I"): 0, ("R2", "Liq", "X_S"): 0, - ("R2", "Liq", "X_BH"): 1, - ("R2", "Liq", "X_BA"): 0, - ("R2", "Liq", "X_P"): 0, - ("R2", "Liq", "S_O"): 0, - ("R2", "Liq", "S_NO"): -(1 - self.Y_H) / (2.86 * self.Y_H), - ("R2", "Liq", "S_NH"): -self.i_xb, - ("R2", "Liq", "S_ND"): 0, - ("R2", "Liq", "X_ND"): 0, - ("R2", "Liq", "S_ALK"): ((1 - self.Y_H) / (2.86 * self.Y_H) - self.i_xb) - * (mw_alk / mw_n), - # R3: Aerobic growth of autotrophs + ("R2", "Liq", "X_H"): 0, + ("R2", "Liq", "X_STO"): self.Y_STO_O2, + ("R2", "Liq", "X_A"): 0, + ("R2", "Liq", "X_TSS"): t2, + # R3: Anoxic storage of S_S ("R3", "Liq", "H2O"): 0, + ("R3", "Liq", "S_O"): 0, ("R3", "Liq", "S_I"): 0, - ("R3", "Liq", "S_S"): 0, + ("R3", "Liq", "S_S"): -1, + ("R3", "Liq", "S_NH4"): y3, + ("R3", "Liq", "S_N2"): -x3, + ("R3", "Liq", "S_NOX"): x3, + ("R3", "Liq", "S_ALK"): z3, ("R3", "Liq", "X_I"): 0, ("R3", "Liq", "X_S"): 0, - ("R3", "Liq", "X_BH"): 0, - ("R3", "Liq", "X_BA"): 1, - ("R3", "Liq", "X_P"): 0, - ("R3", "Liq", "S_O"): -(4.57 - self.Y_A) / self.Y_A, - ("R3", "Liq", "S_NO"): 1 / self.Y_A, - ("R3", "Liq", "S_NH"): -self.i_xb - 1 / self.Y_A, - ("R3", "Liq", "S_ND"): 0, - ("R3", "Liq", "X_ND"): 0, - ("R3", "Liq", "S_ALK"): (-self.i_xb - 2 / (self.Y_A)) * (mw_alk / mw_n), - # R4: Decay of heterotrophs + ("R3", "Liq", "X_H"): 0, + ("R3", "Liq", "X_STO"): self.Y_STO_NOX, + ("R3", "Liq", "X_A"): 0, + ("R3", "Liq", "X_TSS"): t3, + # R4: Aerobic growth of X_H ("R4", "Liq", "H2O"): 0, + ("R4", "Liq", "S_O"): x4, ("R4", "Liq", "S_I"): 0, ("R4", "Liq", "S_S"): 0, + ("R4", "Liq", "S_NH4"): y4, + ("R4", "Liq", "S_N2"): 0, + ("R4", "Liq", "S_NOX"): 0, + ("R4", "Liq", "S_ALK"): z4, ("R4", "Liq", "X_I"): 0, - ("R4", "Liq", "X_S"): 1 - self.f_p, - ("R4", "Liq", "X_BH"): -1, - ("R4", "Liq", "X_BA"): 0, - ("R4", "Liq", "X_P"): self.f_p, - ("R4", "Liq", "S_O"): 0, - ("R4", "Liq", "S_NO"): 0, - ("R4", "Liq", "S_NH"): 0, - ("R4", "Liq", "S_ND"): 0, - ("R4", "Liq", "X_ND"): self.i_xb - self.f_p * self.i_xp, - ("R4", "Liq", "S_ALK"): 0, - # R5: Decay of autotrophs + ("R4", "Liq", "X_S"): 0, + ("R4", "Liq", "X_H"): 1, + ("R4", "Liq", "X_STO"): -1 / self.Y_H_O2, + ("R4", "Liq", "X_A"): 0, + ("R4", "Liq", "X_TSS"): t4, + # R5: Anoxic growth (denitrific.) ("R5", "Liq", "H2O"): 0, + ("R5", "Liq", "S_O"): 0, ("R5", "Liq", "S_I"): 0, ("R5", "Liq", "S_S"): 0, + ("R5", "Liq", "S_NH4"): y4, + ("R5", "Liq", "S_N2"): -x5, + ("R5", "Liq", "S_NOX"): x5, + ("R5", "Liq", "S_ALK"): z5, ("R5", "Liq", "X_I"): 0, - ("R5", "Liq", "X_S"): 1 - self.f_p, - ("R5", "Liq", "X_BH"): 0, - ("R5", "Liq", "X_BA"): -1, - ("R5", "Liq", "X_P"): self.f_p, - ("R5", "Liq", "S_O"): 0, - ("R5", "Liq", "S_NO"): 0, - ("R5", "Liq", "S_NH"): 0, - ("R5", "Liq", "S_ND"): 0, - ("R5", "Liq", "X_ND"): self.i_xb - self.f_p * self.i_xp, - ("R5", "Liq", "S_ALK"): 0, - # R6: Ammonification of soluble organic nitrogen + ("R5", "Liq", "X_S"): 0, + ("R5", "Liq", "X_H"): 1, + ("R5", "Liq", "X_STO"): -1 / self.Y_H_NOX, + ("R5", "Liq", "X_A"): 0, + ("R5", "Liq", "X_TSS"): t5, + # R6: Aerobic endog. respiration ("R6", "Liq", "H2O"): 0, + ("R6", "Liq", "S_O"): x6, ("R6", "Liq", "S_I"): 0, ("R6", "Liq", "S_S"): 0, - ("R6", "Liq", "X_I"): 0, + ("R6", "Liq", "S_NH4"): y6, + ("R6", "Liq", "S_N2"): 0, + ("R6", "Liq", "S_NOX"): 0, + ("R6", "Liq", "S_ALK"): z6, + ("R6", "Liq", "X_I"): self.f_XI, ("R6", "Liq", "X_S"): 0, - ("R6", "Liq", "X_BH"): 0, - ("R6", "Liq", "X_BA"): 0, - ("R6", "Liq", "X_P"): 0, - ("R6", "Liq", "S_O"): 0, - ("R6", "Liq", "S_NO"): 0, - ("R6", "Liq", "S_NH"): 1, - ("R6", "Liq", "S_ND"): -1, - ("R6", "Liq", "X_ND"): 0, - ("R6", "Liq", "S_ALK"): 1 * (mw_alk / mw_n), - # R7: Hydrolysis of entrapped organics + ("R6", "Liq", "X_H"): -1, + ("R6", "Liq", "X_STO"): 0, + ("R6", "Liq", "X_A"): 0, + ("R6", "Liq", "X_TSS"): t6, + # R7: Anoxic endog. respiration ("R7", "Liq", "H2O"): 0, - ("R7", "Liq", "S_I"): 0, - ("R7", "Liq", "S_S"): 1, - ("R7", "Liq", "X_I"): 0, - ("R7", "Liq", "X_S"): -1, - ("R7", "Liq", "X_BH"): 0, - ("R7", "Liq", "X_BA"): 0, - ("R7", "Liq", "X_P"): 0, ("R7", "Liq", "S_O"): 0, - ("R7", "Liq", "S_NO"): 0, - ("R7", "Liq", "S_NH"): 0, - ("R7", "Liq", "S_ND"): 0, - ("R7", "Liq", "X_ND"): 0, - ("R7", "Liq", "S_ALK"): 0, - # R8: Hydrolysis of entrapped organic nitrogen + ("R7", "Liq", "S_I"): 0, + ("R7", "Liq", "S_S"): 0, + ("R7", "Liq", "S_NH4"): y7, + ("R7", "Liq", "S_N2"): -x7, + ("R7", "Liq", "S_NOX"): x7, + ("R7", "Liq", "S_ALK"): z7, + ("R7", "Liq", "X_I"): self.f_XI, + ("R7", "Liq", "X_S"): 0, + ("R7", "Liq", "X_H"): -1, + ("R7", "Liq", "X_STO"): 0, + ("R7", "Liq", "X_A"): 0, + ("R7", "Liq", "X_TSS"): t7, + # R8: Aerobic respiration of X_STO ("R8", "Liq", "H2O"): 0, + ("R8", "Liq", "S_O"): x8, ("R8", "Liq", "S_I"): 0, ("R8", "Liq", "S_S"): 0, + ("R8", "Liq", "S_NH4"): 0, + ("R8", "Liq", "S_N2"): 0, + ("R8", "Liq", "S_NOX"): 0, + ("R8", "Liq", "S_ALK"): 0, ("R8", "Liq", "X_I"): 0, ("R8", "Liq", "X_S"): 0, - ("R8", "Liq", "X_BH"): 0, - ("R8", "Liq", "X_BA"): 0, - ("R8", "Liq", "X_P"): 0, - ("R8", "Liq", "S_O"): 0, - ("R8", "Liq", "S_NO"): 0, - ("R8", "Liq", "S_NH"): 0, - ("R8", "Liq", "S_ND"): 1, - ("R8", "Liq", "X_ND"): -1, - ("R8", "Liq", "S_ALK"): 0, + ("R8", "Liq", "X_H"): 0, + ("R8", "Liq", "X_STO"): -1, + ("R8", "Liq", "X_A"): 0, + ("R8", "Liq", "X_TSS"): t8, + # R9: Anoxic respiration of X_STO + ("R9", "Liq", "H2O"): 0, + ("R9", "Liq", "S_O"): 0, + ("R9", "Liq", "S_I"): 0, + ("R9", "Liq", "S_S"): 0, + ("R9", "Liq", "S_NH4"): 0, + ("R9", "Liq", "S_N2"): -x9, + ("R9", "Liq", "S_NOX"): x9, + ("R9", "Liq", "S_ALK"): z9, + ("R9", "Liq", "X_I"): 0, + ("R9", "Liq", "X_S"): 0, + ("R9", "Liq", "X_H"): 0, + ("R9", "Liq", "X_STO"): -1, + ("R9", "Liq", "X_A"): 0, + ("R9", "Liq", "X_TSS"): t9, + # Autotrophic organisms, nitrifying activity + # R10: Aerobic growth of X_A + ("R10", "Liq", "H2O"): 0, + ("R10", "Liq", "S_O"): x10, + ("R10", "Liq", "S_I"): 0, + ("R10", "Liq", "S_S"): 0, + ("R10", "Liq", "S_NH4"): y10, + ("R10", "Liq", "S_N2"): 0, + ("R10", "Liq", "S_NOX"): 1 / self.Y_A, + ("R10", "Liq", "S_ALK"): z10, + ("R10", "Liq", "X_I"): 0, + ("R10", "Liq", "X_S"): 0, + ("R10", "Liq", "X_H"): 0, + ("R10", "Liq", "X_STO"): 0, + ("R10", "Liq", "X_A"): 1, + ("R10", "Liq", "X_TSS"): t10, + # R11: Aerobic endog. respiration + ("R11", "Liq", "H2O"): 0, + ("R11", "Liq", "S_O"): x11, + ("R11", "Liq", "S_I"): 0, + ("R11", "Liq", "S_S"): 0, + ("R11", "Liq", "S_NH4"): y11, + ("R11", "Liq", "S_N2"): 0, + ("R11", "Liq", "S_NOX"): 0, + ("R11", "Liq", "S_ALK"): z11, + ("R11", "Liq", "X_I"): self.f_XI, + ("R11", "Liq", "X_S"): 0, + ("R11", "Liq", "X_H"): 0, + ("R11", "Liq", "X_STO"): 0, + ("R11", "Liq", "X_A"): -1, + ("R11", "Liq", "X_TSS"): t11, + # R12: Anoxic endog. respiration + ("R12", "Liq", "H2O"): 0, + ("R12", "Liq", "S_O"): 0, + ("R12", "Liq", "S_I"): 0, + ("R12", "Liq", "S_S"): 0, + ("R12", "Liq", "S_NH4"): y12, + ("R12", "Liq", "S_N2"): -x12, + ("R12", "Liq", "S_NOX"): x12, + ("R12", "Liq", "S_ALK"): z12, + ("R12", "Liq", "X_I"): self.f_XI, + ("R12", "Liq", "X_S"): 0, + ("R12", "Liq", "X_H"): 0, + ("R12", "Liq", "X_STO"): 0, + ("R12", "Liq", "X_A"): -1, + ("R12", "Liq", "X_TSS"): t12, } + # Fix all the variables we just created for v in self.component_objects(pyo.Var, descend_into=False): v.fix() @@ -472,111 +643,201 @@ def _rxn_rate(self): def rate_expression_rule(b, r): if r == "R1": - # R1: Aerobic growth of heterotrophs + # R1: Hydrolysis return b.reaction_rate[r] == pyo.units.convert( - b.params.mu_H - * ( - b.conc_mass_comp_ref["S_S"] - / (b.params.K_S + b.conc_mass_comp_ref["S_S"]) + b.params.k_H["20C"] + * (b.conc_mass_comp_ref["X_S"] / b.conc_mass_comp_ref["X_H"]) + / ( + b.params.K_X + + b.conc_mass_comp_ref["X_S"] / b.conc_mass_comp_ref["X_H"] ) + * b.conc_mass_comp_ref["X_H"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + # Heterotrophic organisms, aerobic and denitrifying activity + elif r == "R2": + # R2: Aerobic storage of S_S + return b.reaction_rate[r] == pyo.units.convert( + b.params.k_STO["20C"] * ( b.conc_mass_comp_ref["S_O"] - / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) ) - * b.conc_mass_comp_ref["X_BH"], + * ( + b.conc_mass_comp_ref["S_S"] + / (b.params.K_S + b.conc_mass_comp_ref["S_S"]) + ) + * b.conc_mass_comp_ref["X_H"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R2": - # R2: Anoxic growth of heterotrophs + elif r == "R3": + # R3: Anoxic storage of S_S return b.reaction_rate[r] == pyo.units.convert( - b.params.mu_H + b.params.k_STO["20C"] + * b.params.eta_NOX + * ( + b.params.K_O2 + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) + ) + * ( + b.conc_mass_comp_ref["S_NOX"] + / (b.params.K_NOX + b.conc_mass_comp_ref["S_NOX"]) + ) * ( b.conc_mass_comp_ref["S_S"] / (b.params.K_S + b.conc_mass_comp_ref["S_S"]) ) + * b.conc_mass_comp_ref["X_H"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R4": + # R4: Aerobic growth + return b.reaction_rate[r] == pyo.units.convert( + b.params.mu_H["20C"] + * ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) + ) * ( - b.params.K_OH - / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) + b.conc_mass_comp_ref["S_NH4"] + / (b.params.K_NH4 + b.conc_mass_comp_ref["S_NH4"]) ) * ( - b.conc_mass_comp_ref["S_NO"] - / (b.params.K_NO + b.conc_mass_comp_ref["S_NO"]) + b.conc_mass_comp_ref["S_ALK"] + / (b.params.K_ALK + b.conc_mass_comp_ref["S_ALK"]) ) - * b.params.eta_g - * b.conc_mass_comp_ref["X_BH"], + * (b.conc_mass_comp_ref["X_STO"] / b.conc_mass_comp_ref["X_H"]) + / ( + b.params.K_STO + + b.conc_mass_comp_ref["X_STO"] + / b.conc_mass_comp_ref["X_H"] + ) + * b.conc_mass_comp_ref["X_H"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R3": - # R3: Aerobic growth of autotrophs + elif r == "R5": + # R5: Anoxic growth return b.reaction_rate[r] == pyo.units.convert( - b.params.mu_A + b.params.mu_H["20C"] + * b.params.eta_NOX + * ( + b.params.K_O2 + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) + ) * ( - b.conc_mass_comp_ref["S_NH"] - / (b.params.K_NH + b.conc_mass_comp_ref["S_NH"]) + b.conc_mass_comp_ref["S_NH4"] + / (b.params.K_NH4 + b.conc_mass_comp_ref["S_NH4"]) ) + * ( + b.conc_mass_comp_ref["S_ALK"] + / (b.params.K_ALK + b.conc_mass_comp_ref["S_ALK"]) + ) + * (b.conc_mass_comp_ref["X_STO"] / b.conc_mass_comp_ref["X_H"]) + / ( + b.params.K_STO + + b.conc_mass_comp_ref["X_STO"] + / b.conc_mass_comp_ref["X_H"] + ) + * b.conc_mass_comp_ref["X_H"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R6": + # R6: Aerobic endogenous respiration + return b.reaction_rate[r] == pyo.units.convert( + b.params.b_H_O2["20C"] * ( b.conc_mass_comp_ref["S_O"] - / (b.params.K_OA + b.conc_mass_comp_ref["S_O"]) + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) ) - * b.conc_mass_comp_ref["X_BA"], + * b.conc_mass_comp_ref["X_H"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R4": - # R4: Decay of heterotrophs + elif r == "R7": + # R7: Anoxic endogenous respiration return b.reaction_rate[r] == pyo.units.convert( - b.params.b_H * b.conc_mass_comp_ref["X_BH"], + b.params.b_H_NOX["20C"] + * ( + b.params.K_O2 + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) + ) + * ( + b.conc_mass_comp_ref["S_NOX"] + / (b.params.K_NOX + b.conc_mass_comp_ref["S_NOX"]) + ) + * b.conc_mass_comp_ref["X_H"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R5": - # R5: Decay of autotrophs + elif r == "R8": + # R8: Aerobic respiration of X_STO return b.reaction_rate[r] == pyo.units.convert( - b.params.b_A * b.conc_mass_comp_ref["X_BA"], + b.params.b_STO_O2["20C"] + * ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) + ) + * b.conc_mass_comp_ref["X_STO"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R6": - # R6: Ammonification of soluble organic nitrogen + elif r == "R9": + # R9: Anoxic respiration of X_STO return b.reaction_rate[r] == pyo.units.convert( - b.params.k_a - * b.conc_mass_comp_ref["S_ND"] - * b.conc_mass_comp_ref["X_BH"], + b.params.b_STO_NOX["20C"] + * ( + b.params.K_O2 + / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) + ) + * ( + b.conc_mass_comp_ref["S_NOX"] + / (b.params.K_NOX + b.conc_mass_comp_ref["S_NOX"]) + ) + * b.conc_mass_comp_ref["X_STO"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R7": - # R7: Hydrolysis of entrapped organics + # Autotrophic organisms, nitrifying activity + elif r == "R10": + # R10: Aerobic growth of X_A, nitrification return b.reaction_rate[r] == pyo.units.convert( - b.params.k_h - * b.conc_mass_comp_ref["X_S"] - / ( - b.params.K_X * b.conc_mass_comp_ref["X_BH"] - + b.conc_mass_comp_ref["X_S"] + b.params.mu_A["20C"] + * ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_A_O2 + b.conc_mass_comp_ref["S_O"]) + ) + * ( + b.conc_mass_comp_ref["S_NH4"] + / (b.params.K_A_NH4 + b.conc_mass_comp_ref["S_NH4"]) ) * ( - ( - b.conc_mass_comp_ref["S_O"] - / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) - ) - + b.params.eta_h - * b.params.K_OH - / (b.params.K_OH + b.conc_mass_comp_ref["S_O"]) - * ( - b.conc_mass_comp_ref["S_NO"] - / (b.params.K_NO + b.conc_mass_comp_ref["S_NO"]) - ) + b.conc_mass_comp_ref["S_ALK"] + / (b.params.K_A_ALK + b.conc_mass_comp_ref["S_ALK"]) ) - * b.conc_mass_comp_ref["X_BH"], + * b.conc_mass_comp_ref["X_A"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) - elif r == "R8": - # R8: Hydrolysis of entrapped organic nitrogen - return b.reaction_rate[r] == ( - b.reaction_rate["R7"] + elif r == "R11": + # R11: Aerobic endogenous respiration + return b.reaction_rate[r] == pyo.units.convert( + b.params.b_A_O2["20C"] + * ( + b.conc_mass_comp_ref["S_O"] + / (b.params.K_A_O2 + b.conc_mass_comp_ref["S_O"]) + ) + * b.conc_mass_comp_ref["X_A"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, + ) + elif r == "R12": + # R12: Anoxic endogenous respiration + return b.reaction_rate[r] == pyo.units.convert( + b.params.b_A_NOX["20C"] * ( - b.conc_mass_comp_ref["X_ND"] - / ( - b.conc_mass_comp_ref["X_S"] - + 1e-10 * pyo.units.kg / pyo.units.m**3 - ) + b.params.K_A_O2 + / (b.params.K_A_O2 + b.conc_mass_comp_ref["S_O"]) ) + * ( + b.conc_mass_comp_ref["S_NOX"] + / (b.params.K_A_NOX + b.conc_mass_comp_ref["S_NOX"]) + ) + * b.conc_mass_comp_ref["X_A"], + to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, ) else: raise BurntToast() @@ -598,6 +859,6 @@ def get_reaction_rate_basis(self): def calculate_scaling_factors(self): super().calculate_scaling_factors() - iscale.constraint_scaling_transform(self.rate_expression["R5"], 1e3) - iscale.constraint_scaling_transform(self.rate_expression["R3"], 1e3) - iscale.constraint_scaling_transform(self.rate_expression["R4"], 1e3) + # iscale.constraint_scaling_transform(self.rate_expression["R5"], 1e3) + # iscale.constraint_scaling_transform(self.rate_expression["R3"], 1e3) + # iscale.constraint_scaling_transform(self.rate_expression["R4"], 1e3) From 3df82248b0965cec777e3ea9362d917722b820e9 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Thu, 12 Jun 2025 10:27:26 -0400 Subject: [PATCH 08/21] add test file --- .../activated_sludge/asm3_reactions.py | 2 +- .../tests/test_asm3_reaction.py | 381 ++++++++++++++++++ 2 files changed, 382 insertions(+), 1 deletion(-) create mode 100644 watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py index 2c536c5af8..2c632357d6 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py @@ -367,7 +367,7 @@ def build(self): ("R1", "Liq", "X_H"): 0, ("R1", "Liq", "X_STO"): 0, ("R1", "Liq", "X_A"): 0, - ("R1", "Liq", "X_TSS"): -self.i_SSXS, + ("R1", "Liq", "X_TSS"): t1, # Heterotrophic organisms, aerobic and denitrifying activity # R2: Aerobic storage of S_S ("R2", "Liq", "H2O"): 0, diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py new file mode 100644 index 0000000000..2c1d35a26d --- /dev/null +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py @@ -0,0 +1,381 @@ +################################################################################# +# 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/" +################################################################################# +""" +Tests for ASM3 reaction package. + +Authors: Chenyu Wang +""" +import pytest + +from pyomo.environ import ( + check_optimal_termination, + ConcreteModel, + Constraint, + Suffix, + units, + value, + Var, +) +from pyomo.util.check_units import assert_units_consistent + +from idaes.core import FlowsheetBlock +from idaes.models.unit_models import CSTR +from idaes.core import MaterialFlowBasis +from watertap.core.solvers import get_solver +from idaes.core.util.model_statistics import degrees_of_freedom + +from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( + ASM3ParameterBlock, +) +from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( + ASM3ReactionParameterBlock, + ASM3ReactionBlock, + ASM3ReactionScaler, +) + +# ----------------------------------------------------------------------------- +# Get default solver for testing +solver = get_solver() + + +class TestParamBlock(object): + @pytest.fixture(scope="class") + def model(self): + model = ConcreteModel() + model.pparams = ASM3ParameterBlock() + model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) + + return model + + +# @pytest.mark.unit +# def test_config(self, model): +# assert len(model.rparams.config) == 2 +# +# @pytest.mark.unit +# def test_build(self, model): +# assert model.rparams.reaction_block_class is ASM3ReactionBlock +# +# assert len(model.rparams.rate_reaction_idx) == 8 +# for i in model.rparams.rate_reaction_idx: +# assert i in ["R1", "R2", "R3", "R4", "R5", "R6", "R7", "R8"] +# +# assert len(model.rparams.rate_reaction_stoichiometry) == 8 * 14 +# for i in model.rparams.rate_reaction_stoichiometry: +# assert i[0] in ["R1", "R2", "R3", "R4", "R5", "R6", "R7", "R8"] +# assert i[1] == "Liq" +# assert i[2] in [ +# "H2O", +# "S_I", +# "S_S", +# "X_I", +# "X_S", +# "X_BH", +# "X_BA", +# "X_P", +# "S_O", +# "S_NO", +# "S_NH", +# "S_ND", +# "X_ND", +# "S_ALK", +# ] +# +# assert isinstance(model.rparams.Y_A, Var) +# assert value(model.rparams.Y_A) == 0.24 +# assert isinstance(model.rparams.Y_H, Var) +# assert value(model.rparams.Y_H) == 0.67 +# assert isinstance(model.rparams.f_p, Var) +# assert value(model.rparams.f_p) == 0.08 +# assert isinstance(model.rparams.i_xb, Var) +# assert value(model.rparams.i_xb) == 0.08 +# assert isinstance(model.rparams.i_xp, Var) +# assert value(model.rparams.i_xp) == 0.06 +# +# assert isinstance(model.rparams.mu_A, Var) +# assert value(model.rparams.mu_A) == 0.5 +# assert isinstance(model.rparams.mu_H, Var) +# assert value(model.rparams.mu_H) == 4 +# assert isinstance(model.rparams.K_S, Var) +# assert value(model.rparams.K_S) == 10e-3 +# assert isinstance(model.rparams.K_OH, Var) +# assert value(model.rparams.K_OH) == 0.2e-3 +# assert isinstance(model.rparams.K_OA, Var) +# assert value(model.rparams.K_OA) == 0.4e-3 +# assert isinstance(model.rparams.K_NO, Var) +# assert value(model.rparams.K_NO) == 0.5e-3 +# assert isinstance(model.rparams.b_H, Var) +# assert value(model.rparams.b_H) == 0.3 +# assert isinstance(model.rparams.b_A, Var) +# assert value(model.rparams.b_A) == 0.05 +# assert isinstance(model.rparams.eta_g, Var) +# assert value(model.rparams.eta_g) == 0.8 +# assert isinstance(model.rparams.eta_h, Var) +# assert value(model.rparams.eta_h) == 0.8 +# assert isinstance(model.rparams.k_h, Var) +# assert value(model.rparams.k_h) == 3 +# assert isinstance(model.rparams.K_X, Var) +# assert value(model.rparams.K_X) == 0.1 +# assert isinstance(model.rparams.K_NH, Var) +# assert value(model.rparams.K_NH) == 1e-3 +# assert isinstance(model.rparams.k_a, Var) +# assert value(model.rparams.k_a) == 50 +# +# +# class TestReactionBlock(object): +# @pytest.fixture(scope="class") +# def model(self): +# model = ConcreteModel() +# model.pparams = ASM3ParameterBlock() +# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) +# +# model.props = model.pparams.build_state_block([1]) +# +# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) +# +# return model +# +# @pytest.mark.unit +# def test_build(self, model): +# assert model.rxns[1].conc_mass_comp_ref is model.props[1].conc_mass_comp +# +# @pytest.mark.unit +# def test_rxn_rate(self, model): +# assert isinstance(model.rxns[1].reaction_rate, Var) +# assert len(model.rxns[1].reaction_rate) == 8 +# assert isinstance(model.rxns[1].rate_expression, Constraint) +# assert len(model.rxns[1].rate_expression) == 8 +# +# @pytest.mark.unit +# def test_get_reaction_rate_basis(self, model): +# assert model.rxns[1].get_reaction_rate_basis() == MaterialFlowBasis.mass +# +# @pytest.mark.component +# def test_initialize(self, model): +# assert model.rxns.initialize() is None +# +# @pytest.mark.unit +# def check_units(self, model): +# assert_units_consistent(model) +# +# +# class TestASM3ReactionScaler(object): +# @pytest.mark.unit +# def test_variable_scaling_routine(self): +# model = ConcreteModel() +# model.pparams = ASM3ParameterBlock() +# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) +# +# model.props = model.pparams.build_state_block([1]) +# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) +# +# # Trigger build of reaction properties +# model.rxns[1].reaction_rate +# +# scaler = model.rxns[1].default_scaler() +# assert isinstance(scaler, ASM3ReactionScaler) +# +# scaler.variable_scaling_routine(model.rxns[1]) +# +# assert isinstance(model.rxns[1].scaling_factor, Suffix) +# +# sfx = model.rxns[1].scaling_factor +# assert len(sfx) == 8 +# assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8) +# +# @pytest.mark.unit +# def test_constraint_scaling_routine(self): +# model = ConcreteModel() +# model.pparams = ASM3ParameterBlock() +# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) +# +# model.props = model.pparams.build_state_block([1]) +# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) +# +# # Trigger build of reaction properties +# model.rxns[1].reaction_rate +# +# scaler = model.rxns[1].default_scaler() +# assert isinstance(scaler, ASM3ReactionScaler) +# +# scaler.constraint_scaling_routine(model.rxns[1]) +# +# assert isinstance(model.rxns[1].scaling_factor, Suffix) +# +# sfx = model.rxns[1].scaling_factor +# assert len(sfx) == 8 +# assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx( +# 2.380752e5, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx( +# 1.49540985e8, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx( +# 1.75226112e6, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx( +# 2.88e6, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx( +# 1.728e7, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx( +# 1.728e5, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx( +# 3.174336e5, rel=1e-8 +# ) +# assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1, rel=1e-8) +# +# @pytest.mark.unit +# def test_scale_model(self): +# model = ConcreteModel() +# model.pparams = ASM3ParameterBlock() +# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) +# +# model.props = model.pparams.build_state_block([1]) +# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) +# +# # Trigger build of reaction properties +# model.rxns[1].reaction_rate +# +# scaler = model.rxns[1].default_scaler() +# assert isinstance(scaler, ASM3ReactionScaler) +# +# scaler.scale_model(model.rxns[1]) +# +# assert isinstance(model.rxns[1].scaling_factor, Suffix) +# +# sfx = model.rxns[1].scaling_factor +# assert len(sfx) == 16 +# assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8) +# +# assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx(1e2, rel=1e-8) +# assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1e2, rel=1e-8) +# +# +# class TestReactor: +# @pytest.fixture(scope="class") +# def model(self): +# m = ConcreteModel() +# +# m.fs = FlowsheetBlock(dynamic=False) +# +# m.fs.props = ASM3ParameterBlock() +# m.fs.rxn_props = ASM3ReactionParameterBlock(property_package=m.fs.props) +# +# m.fs.R1 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) +# +# # Feed conditions based on manual mass balance of inlet and recycle streams +# m.fs.R1.inlet.flow_vol.fix(92230 * units.m**3 / units.day) +# m.fs.R1.inlet.temperature.fix(298.15 * units.K) +# m.fs.R1.inlet.pressure.fix(1 * units.atm) +# m.fs.R1.inlet.conc_mass_comp[0, "S_I"].fix(30 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "S_S"].fix(14.6112 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "X_I"].fix(1149 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "X_S"].fix(89.324 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "X_BH"].fix(2542.03 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "X_BA"].fix(148.6 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "X_P"].fix(448 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "S_O"].fix(0.3928 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "S_NO"].fix(8.32 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "S_NH"].fix(7.696 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "S_ND"].fix(1.9404 * units.g / units.m**3) +# m.fs.R1.inlet.conc_mass_comp[0, "X_ND"].fix(5.616 * units.g / units.m**3) +# m.fs.R1.inlet.alkalinity.fix(4.704 * units.mol / units.m**3) +# +# m.fs.R1.volume.fix(1000 * units.m**3) +# +# return m +# +# @pytest.mark.unit +# def test_dof(self, model): +# assert degrees_of_freedom(model) == 0 +# +# @pytest.mark.unit +# def test_unit_consistency(self, model): +# assert_units_consistent(model) == 0 +# +# @pytest.mark.component +# def test_solve(self, model): +# model.fs.R1.initialize() +# +# solver = get_solver() +# results = solver.solve(model, tee=True) +# assert check_optimal_termination(results) +# +# @pytest.mark.component +# def test_solution(self, model): +# assert value(model.fs.R1.outlet.flow_vol[0]) == pytest.approx(1.0675, rel=1e-4) +# assert value(model.fs.R1.outlet.temperature[0]) == pytest.approx( +# 298.15, rel=1e-4 +# ) +# assert value(model.fs.R1.outlet.pressure[0]) == pytest.approx(101325, rel=1e-4) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_I"]) == pytest.approx( +# 30e-3, rel=1e-5 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_S"]) == pytest.approx( +# 2.81e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_I"]) == pytest.approx( +# 1149e-3, rel=1e-3 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_S"]) == pytest.approx( +# 82.1e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_BH"]) == pytest.approx( +# 2552e-3, rel=1e-3 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_BA"]) == pytest.approx( +# 149e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_P"]) == pytest.approx( +# 449e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_O"]) == pytest.approx( +# 4.3e-6, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NO"]) == pytest.approx( +# 5.36e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NH"]) == pytest.approx( +# 7.92e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_ND"]) == pytest.approx( +# 1.22e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_ND"]) == pytest.approx( +# 5.29e-3, rel=1e-2 +# ) +# assert value(model.fs.R1.outlet.alkalinity[0]) == pytest.approx( +# 4.93e-3, rel=1e-2 +# ) From 03683301899abb6a726c83dce327cfc9913d7424 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Thu, 12 Jun 2025 16:56:03 -0400 Subject: [PATCH 09/21] add test file --- .../tests/test_asm3_reaction.py | 4 +- .../tests/test_asm3_thermo.py | 404 ++++++++++++++++++ 2 files changed, 406 insertions(+), 2 deletions(-) create mode 100644 watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py index 2c1d35a26d..e66189a150 100644 --- a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py @@ -33,10 +33,10 @@ from watertap.core.solvers import get_solver from idaes.core.util.model_statistics import degrees_of_freedom -from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( +from watertap.property_models.unit_specific.activated_sludge.asm3_properties import ( ASM3ParameterBlock, ) -from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( +from watertap.property_models.unit_specific.activated_sludge.asm3_reactions import ( ASM3ReactionParameterBlock, ASM3ReactionBlock, ASM3ReactionScaler, diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py new file mode 100644 index 0000000000..75477d8f3a --- /dev/null +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py @@ -0,0 +1,404 @@ +################################################################################# +# 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/" +################################################################################# +""" +Tests for ASM3 thermo property package. +Authors: Chenyu Wang +""" + +import pytest +from pyomo.environ import ConcreteModel, Expression, Param, units, value, Var, Suffix +from pyomo.util.check_units import assert_units_consistent +from idaes.core import MaterialBalanceType, EnergyBalanceType, MaterialFlowBasis + +from watertap.property_models.unit_specific.activated_sludge.asm3_properties import ( + ASM3ParameterBlock, + ASM3StateBlock, + ASM3PropertiesScaler, +) +from idaes.core.util.model_statistics import ( + fixed_variables_set, + activated_constraints_set, +) + +from watertap.core.solvers import get_solver + + +# ----------------------------------------------------------------------------- +# Get default solver for testing +solver = get_solver() + + +class TestParamBlock(object): + @pytest.fixture(scope="class") + def model(self): + model = ConcreteModel() + model.params = ASM3ParameterBlock() + + return model + + @pytest.mark.unit + def test_build(self, model): + assert model.params.state_block_class is ASM3StateBlock + + assert len(model.params.phase_list) == 1 + for i in model.params.phase_list: + assert i == "Liq" + + assert len(model.params.component_list) == 14 + for i in model.params.component_list: + assert i in [ + "H2O", + "S_O", + "S_I", + "S_S", + "X_S", + "S_NH4", + "S_N2", + "S_NOX", + "S_ALK", + "X_I", + "X_S", + "X_H", + "X_STO", + "X_A", + "X_TSS", + ] + # + assert isinstance(model.params.cp_mass, Param) + assert value(model.params.cp_mass) == 4182 + + assert isinstance(model.params.dens_mass, Param) + assert value(model.params.dens_mass) == 997 + + assert isinstance(model.params.pressure_ref, Param) + assert value(model.params.pressure_ref) == 101325 + + assert isinstance(model.params.temperature_ref, Param) + assert value(model.params.temperature_ref) == 298.15 + + assert len(model.params.particulate_component_set) == 6 + assert len(model.params.non_particulate_component_set) == 8 + for i in model.params.particulate_component_set: + assert i in [ + "X_I", + "X_S", + "X_H", + "X_STO", + "X_A", + "X_TSS", + ] + + for i in model.params.non_particulate_component_set: + assert i in [ + "H2O", + "S_O", + "S_I", + "S_S", + "X_S", + "S_NH4", + "S_N2", + "S_NOX", + "S_ALK", + ] + + +class TestStateBlock(object): + @pytest.fixture(scope="class") + def model(self): + model = ConcreteModel() + model.params = ASM3ParameterBlock() + + model.props = model.params.build_state_block([1]) + + return model + + @pytest.mark.unit + def test_build(self, model): + assert model.props[1].default_scaler is ASM3PropertiesScaler + + assert isinstance(model.props[1].flow_vol, Var) + assert value(model.props[1].flow_vol) == 1 + + assert isinstance(model.props[1].pressure, Var) + assert value(model.props[1].pressure) == 101325 + + assert isinstance(model.props[1].temperature, Var) + assert value(model.props[1].temperature) == 298.15 + + assert isinstance(model.props[1].alkalinity, Var) + assert value(model.props[1].alkalinity) == 1 + + assert isinstance(model.props[1].conc_mass_comp, Var) + # H2O should not appear in conc_mass_comp + assert len(model.props[1].conc_mass_comp) == 12 + for i in model.props[1].conc_mass_comp: + assert i in [ + "S_O", + "S_I", + "S_S", + "X_S", + "S_NH4", + "S_N2", + "S_NOX", + "S_ALK", + "X_I", + "X_S", + "X_H", + "X_STO", + "X_A", + "X_TSS", + ] + assert value(model.props[1].conc_mass_comp[i]) == 0.1 + + assert isinstance(model.props[1].params.f_SI, Var) + assert value(model.props[1].params.f_SI) == 0 + assert isinstance(model.props[1].params.f_XI, Var) + assert value(model.props[1].params.f_XI) == 0.2 + assert isinstance(model.props[1].params.i_NSI, Var) + assert value(model.props[1].params.i_NSI) == 0.01 + assert isinstance(model.props[1].params.i_NSS, Var) + assert value(model.props[1].params.i_NSS) == 0.03 + assert isinstance(model.props[1].params.i_NXI, Var) + assert value(model.props[1].params.i_NXI) == 0.02 + assert isinstance(model.props[1].params.i_NXS, Var) + assert value(model.props[1].params.i_NXS) == 0.04 + assert isinstance(model.props[1].params.i_NBM, Var) + assert value(model.props[1].params.i_NBM) == 0.07 + assert isinstance(model.props[1].params.i_SSXI, Var) + assert value(model.props[1].params.i_SSXI) == 0.75 + assert isinstance(model.props[1].params.i_SSXS, Var) + assert value(model.props[1].params.i_SSXS) == 0.75 + assert isinstance(model.props[1].params.i_SSBM, Var) + assert value(model.props[1].params.i_SSBM) == 0.90 + assert isinstance(model.props[1].params.i_SSSTO, Var) + assert value(model.props[1].params.i_SSSTO) == 0.60 + + assert isinstance(model.props[1].material_flow_expression, Expression) + for j in model.params.component_list: + if j == "H2O": + assert str(model.props[1].material_flow_expression[j].expr) == str( + model.props[1].flow_vol * model.props[1].params.dens_mass + ) + elif j == "S_ALK": + assert str(model.props[1].material_flow_expression[j].expr) == str( + model.props[1].flow_vol + * model.props[1].alkalinity + * (12 * units.kg / units.kmol) + ) + else: + assert str(model.props[1].material_flow_expression[j].expr) == str( + model.props[1].flow_vol * model.props[1].conc_mass_comp[j] + ) + + assert isinstance(model.props[1].enthalpy_flow_expression, Expression) + assert str(model.props[1].enthalpy_flow_expression.expr) == str( + model.props[1].flow_vol + * model.props[1].params.dens_mass + * model.props[1].params.cp_mass + * (model.props[1].temperature - model.props[1].params.temperature_ref) + ) + + assert isinstance(model.props[1].material_density_expression, Expression) + for j in model.params.component_list: + if j == "H2O": + assert model.props[1].material_density_expression[j].expr is ( + model.props[1].params.dens_mass + ) + elif j == "S_ALK": + assert str(model.props[1].material_density_expression[j].expr) == str( + model.props[1].alkalinity * (12 * units.kg / units.kmol) + ) + else: + assert str(model.props[1].material_density_expression[j].expr) == str( + model.props[1].conc_mass_comp[j] + ) + + assert isinstance(model.props[1].energy_density_expression, Expression) + assert str(model.props[1].energy_density_expression.expr) == str( + model.props[1].params.dens_mass + * model.props[1].params.cp_mass + * (model.props[1].temperature - model.props[1].params.temperature_ref) + ) + + @pytest.mark.unit + def test_get_material_flow_terms(self, model): + for p in model.params.phase_list: + for j in model.params.component_list: + assert model.props[1].get_material_flow_terms(p, j) is ( + model.props[1].material_flow_expression[j] + ) + + @pytest.mark.unit + def test_get_enthalpy_flow_terms(self, model): + for p in model.params.phase_list: + assert model.props[1].get_enthalpy_flow_terms(p) is ( + model.props[1].enthalpy_flow_expression + ) + + @pytest.mark.unit + def test_get_material_density_terms(self, model): + for p in model.params.phase_list: + for j in model.params.component_list: + assert model.props[1].get_material_density_terms(p, j) is ( + model.props[1].material_density_expression[j] + ) + + @pytest.mark.unit + def test_get_energy_density_terms(self, model): + for p in model.params.phase_list: + assert model.props[1].get_energy_density_terms(p) is ( + model.props[1].energy_density_expression + ) + + @pytest.mark.unit + def test_default_material_balance_type(self, model): + assert ( + model.props[1].default_material_balance_type() + == MaterialBalanceType.componentPhase + ) + + @pytest.mark.unit + def test_default_energy_balance_type(self, model): + assert ( + model.props[1].default_energy_balance_type() + == EnergyBalanceType.enthalpyTotal + ) + + @pytest.mark.unit + def test_get_material_flow_basis(self, model): + assert model.props[1].get_material_flow_basis() == MaterialFlowBasis.mass + + @pytest.mark.unit + def test_define_state_vars(self, model): + sv = model.props[1].define_state_vars() + + assert len(sv) == 5 + for i in sv: + assert i in [ + "flow_vol", + "alkalinity", + "conc_mass_comp", + "temperature", + "pressure", + ] + + @pytest.mark.unit + def test_define_port_members(self, model): + sv = model.props[1].define_port_members() + + assert len(sv) == 5 + for i in sv: + assert i in [ + "flow_vol", + "alkalinity", + "conc_mass_comp", + "temperature", + "pressure", + ] + + @pytest.mark.unit + def test_define_display_vars(self, model): + sv = model.props[1].define_display_vars() + + assert len(sv) == 5 + for i in sv: + assert i in [ + "Volumetric Flowrate", + "Molar Alkalinity", + "Mass Concentration", + "Temperature", + "Pressure", + ] + + @pytest.mark.component + def test_initialize(self, model): + orig_fixed_vars = fixed_variables_set(model) + orig_act_consts = activated_constraints_set(model) + + model.props.initialize(hold_state=False) + + fin_fixed_vars = fixed_variables_set(model) + fin_act_consts = activated_constraints_set(model) + + assert len(fin_act_consts) == len(orig_act_consts) + assert len(fin_fixed_vars) == len(orig_fixed_vars) + + for c in fin_act_consts: + assert c in orig_act_consts + for v in fin_fixed_vars: + assert v in orig_fixed_vars + + @pytest.mark.unit + def check_units(self, model): + assert_units_consistent(model) + + +# @pytest.mark.unit +# def test_expressions(self, model): +# assert value(model.props[1].TSS) == 0.375 +# assert value(model.props[1].COD) == pytest.approx(0.7, rel=1e-3) +# assert value(model.props[1].BOD5["effluent"]) == 0.096 +# assert value(model.props[1].BOD5["raw"]) == 0.096 * 0.65 / 0.25 +# assert value(model.props[1].TKN) == pytest.approx(0.328, rel=1e-3) +# assert value(model.props[1].Total_N) == pytest.approx(0.428, rel=1e-3) + + +class TestASM3PropertiesScaler: + @pytest.mark.unit + def test_variable_scaling_routine(self): + model = ConcreteModel() + model.params = ASM3ParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=False) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, ASM3PropertiesScaler) + + scaler.variable_scaling_routine(model.props[1]) + + sfx = model.props[1].scaling_factor + assert len(sfx) == 3 + assert sfx[model.props[1].flow_vol] == pytest.approx(1e1, rel=1e-8) + assert sfx[model.props[1].pressure] == pytest.approx(1e-6, rel=1e-8) + assert sfx[model.props[1].temperature] == pytest.approx(1e-1, rel=1e-8) + + @pytest.mark.unit + def test_constraint_scaling_routine(self): + model = ConcreteModel() + model.params = ASM3ParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=False) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, ASM3PropertiesScaler) + + scaler.constraint_scaling_routine(model.props[1]) + + @pytest.mark.unit + def test_scale_model(self): + model = ConcreteModel() + model.params = ASM3ParameterBlock() + + model.props = model.params.build_state_block([1], defined_state=False) + + scaler = model.props[1].default_scaler() + assert isinstance(scaler, ASM3PropertiesScaler) + + scaler.scale_model(model.props[1]) + + assert isinstance(model.props[1].scaling_factor, Suffix) + + sfx = model.props[1].scaling_factor + assert len(sfx) == 3 + assert sfx[model.props[1].flow_vol] == pytest.approx(1e1, rel=1e-8) + assert sfx[model.props[1].pressure] == pytest.approx(1e-6, rel=1e-8) + assert sfx[model.props[1].temperature] == pytest.approx(1e-1, rel=1e-8) From bb34640ed4fb73a89b5a856d8025f1e45f18e048 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Mon, 16 Jun 2025 09:53:46 -0400 Subject: [PATCH 10/21] Reaction and properties test pass --- .../activated_sludge/asm3_properties.py | 24 +- .../activated_sludge/asm3_reactions.py | 62 +- .../tests/test_asm3_reaction.py | 742 ++++++++++-------- 3 files changed, 463 insertions(+), 365 deletions(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py index 79521215ac..c1867d9fee 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py @@ -244,16 +244,16 @@ def define_metadata(cls, obj): "conc_mass_comp": {"method": None}, } ) - # obj.define_custom_properties( - # { - # "alkalinity": {"method": None}, - # "TSS": {"method": "_TSS"}, - # "BOD5": {"method": "_BOD5"}, - # "TKN": {"method": "_TKN"}, - # "Total_N": {"method": "_Total_N"}, - # "COD": {"method": "_COD"}, - # } - # ) + obj.define_custom_properties( + { + "alkalinity": {"method": None}, + # "TSS": {"method": "_TSS"}, + # "BOD5": {"method": "_BOD5"}, + # "TKN": {"method": "_TKN"}, + # "Total_N": {"method": "_Total_N"}, + # "COD": {"method": "_COD"}, + } + ) obj.add_default_units( { "time": pyo.units.s, @@ -450,11 +450,11 @@ def material_flow_expression(self, j): if j == "H2O": return self.flow_vol * self.params.dens_mass elif j == "S_ALK": - # Convert moles of alkalinity to mass of C assuming all is HCO3- + # Convert moles of alkalinity to mass assuming all is HCO3- return ( self.flow_vol * self.alkalinity - * (12 * pyo.units.kg / pyo.units.kmol) + * (61 * pyo.units.kg / pyo.units.kmol) ) else: return self.flow_vol * self.conc_mass_comp[j] diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py index 2c632357d6..2217dde9f5 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py @@ -174,22 +174,22 @@ def build(self): doc="Anoxic reduction factor", ) self.K_O2 = pyo.Var( - initialize=0.2, - units=pyo.units.g / pyo.units.m**3, + initialize=0.2e-3, + units=pyo.units.kg / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Saturation constant for S_NO2 (g-O2 / m3)", + doc="Saturation constant for S_NO2 (kg-O2 / m3)", ) self.K_NOX = pyo.Var( - initialize=0.5, - units=pyo.units.g / pyo.units.m**3, + initialize=0.5e-3, + units=pyo.units.kg / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Saturation constant for S_NOX (g-N-NO3 / m3)", + doc="Saturation constant for S_NOX (kg-N-NO3 / m3)", ) self.K_S = pyo.Var( - initialize=2, - units=pyo.units.g / pyo.units.m**3, + initialize=2e-3, + units=pyo.units.kg / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Saturation constant for for substrate S_S (g-COD-S_S / m3)", + doc="Saturation constant for substrate S_S (kg-COD-S_S / m3)", ) self.K_STO = pyo.Var( initialize=1, @@ -206,16 +206,16 @@ def build(self): doc="Heterotrophic max. growth rate of X_H (day^-1)", ) self.K_NH4 = pyo.Var( - initialize=0.01, - units=pyo.units.g / pyo.units.m**3, + initialize=0.01e-3, + units=pyo.units.kg / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Saturation constant for ammonium, S_NH4 (g-N / m3)", + doc="Saturation constant for ammonium, S_NH4 (kg-N / m3)", ) self.K_ALK = pyo.Var( - initialize=0.1, - units=pyo.units.mol / pyo.units.m**3, + initialize=0.1e-3, + units=pyo.units.kmol / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Saturation constant for alkalinity for X_H (mol-HCO3- / m3)", + doc="Saturation constant for alkalinity for X_H (kmol-HCO3- / m3)", ) b_H_O2_dict = {"10C": 0.1, "20C": 0.2} self.b_H_O2 = pyo.Var( @@ -260,22 +260,22 @@ def build(self): doc="Autotrophic max. growth rate of X_A (day^-1)", ) self.K_A_NH4 = pyo.Var( - initialize=1, - units=pyo.units.g / pyo.units.m**3, + initialize=1e-3, + units=pyo.units.kg / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Ammonium substrate saturation for X_A (g-N / m3)", + doc="Ammonium substrate saturation for X_A (kg-N / m3)", ) self.K_A_O2 = pyo.Var( - initialize=0.5, - units=pyo.units.g / pyo.units.m**3, + initialize=0.5e-3, + units=pyo.units.kg / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Oxygen saturation for nitrifiers (g-O2 / m3)", + doc="Oxygen saturation for nitrifiers (kg-O2 / m3)", ) self.K_A_ALK = pyo.Var( - initialize=0.5, - units=pyo.units.mol / pyo.units.m**3, + initialize=0.5e-3, + units=pyo.units.kmol / pyo.units.m**3, domain=pyo.PositiveReals, - doc="Bicarbonate saturation for nitrifiers (mol-HCO3- / m3)", + doc="Bicarbonate saturation for nitrifiers (kmol-HCO3- / m3)", ) b_A_O2_dict = {"10C": 0.05, "20C": 0.15} self.b_A_O2 = pyo.Var( @@ -703,8 +703,8 @@ def rate_expression_rule(b, r): / (b.params.K_NH4 + b.conc_mass_comp_ref["S_NH4"]) ) * ( - b.conc_mass_comp_ref["S_ALK"] - / (b.params.K_ALK + b.conc_mass_comp_ref["S_ALK"]) + b.state_ref.alkalinity + / (b.params.K_ALK + b.state_ref.alkalinity) ) * (b.conc_mass_comp_ref["X_STO"] / b.conc_mass_comp_ref["X_H"]) / ( @@ -729,8 +729,8 @@ def rate_expression_rule(b, r): / (b.params.K_NH4 + b.conc_mass_comp_ref["S_NH4"]) ) * ( - b.conc_mass_comp_ref["S_ALK"] - / (b.params.K_ALK + b.conc_mass_comp_ref["S_ALK"]) + b.state_ref.alkalinity + / (b.params.K_ALK + b.state_ref.alkalinity) ) * (b.conc_mass_comp_ref["X_STO"] / b.conc_mass_comp_ref["X_H"]) / ( @@ -807,8 +807,8 @@ def rate_expression_rule(b, r): / (b.params.K_A_NH4 + b.conc_mass_comp_ref["S_NH4"]) ) * ( - b.conc_mass_comp_ref["S_ALK"] - / (b.params.K_A_ALK + b.conc_mass_comp_ref["S_ALK"]) + b.state_ref.alkalinity + / (b.params.K_A_ALK + b.state_ref.alkalinity) ) * b.conc_mass_comp_ref["X_A"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, @@ -834,7 +834,7 @@ def rate_expression_rule(b, r): ) * ( b.conc_mass_comp_ref["S_NOX"] - / (b.params.K_A_NOX + b.conc_mass_comp_ref["S_NOX"]) + / (b.params.K_NOX + b.conc_mass_comp_ref["S_NOX"]) ) * b.conc_mass_comp_ref["X_A"], to_units=pyo.units.kg / pyo.units.m**3 / pyo.units.s, diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py index e66189a150..bf1850e08c 100644 --- a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py @@ -56,326 +56,424 @@ def model(self): return model + @pytest.mark.unit + def test_config(self, model): + assert len(model.rparams.config) == 2 -# @pytest.mark.unit -# def test_config(self, model): -# assert len(model.rparams.config) == 2 -# -# @pytest.mark.unit -# def test_build(self, model): -# assert model.rparams.reaction_block_class is ASM3ReactionBlock -# -# assert len(model.rparams.rate_reaction_idx) == 8 -# for i in model.rparams.rate_reaction_idx: -# assert i in ["R1", "R2", "R3", "R4", "R5", "R6", "R7", "R8"] -# -# assert len(model.rparams.rate_reaction_stoichiometry) == 8 * 14 -# for i in model.rparams.rate_reaction_stoichiometry: -# assert i[0] in ["R1", "R2", "R3", "R4", "R5", "R6", "R7", "R8"] -# assert i[1] == "Liq" -# assert i[2] in [ -# "H2O", -# "S_I", -# "S_S", -# "X_I", -# "X_S", -# "X_BH", -# "X_BA", -# "X_P", -# "S_O", -# "S_NO", -# "S_NH", -# "S_ND", -# "X_ND", -# "S_ALK", -# ] -# -# assert isinstance(model.rparams.Y_A, Var) -# assert value(model.rparams.Y_A) == 0.24 -# assert isinstance(model.rparams.Y_H, Var) -# assert value(model.rparams.Y_H) == 0.67 -# assert isinstance(model.rparams.f_p, Var) -# assert value(model.rparams.f_p) == 0.08 -# assert isinstance(model.rparams.i_xb, Var) -# assert value(model.rparams.i_xb) == 0.08 -# assert isinstance(model.rparams.i_xp, Var) -# assert value(model.rparams.i_xp) == 0.06 -# -# assert isinstance(model.rparams.mu_A, Var) -# assert value(model.rparams.mu_A) == 0.5 -# assert isinstance(model.rparams.mu_H, Var) -# assert value(model.rparams.mu_H) == 4 -# assert isinstance(model.rparams.K_S, Var) -# assert value(model.rparams.K_S) == 10e-3 -# assert isinstance(model.rparams.K_OH, Var) -# assert value(model.rparams.K_OH) == 0.2e-3 -# assert isinstance(model.rparams.K_OA, Var) -# assert value(model.rparams.K_OA) == 0.4e-3 -# assert isinstance(model.rparams.K_NO, Var) -# assert value(model.rparams.K_NO) == 0.5e-3 -# assert isinstance(model.rparams.b_H, Var) -# assert value(model.rparams.b_H) == 0.3 -# assert isinstance(model.rparams.b_A, Var) -# assert value(model.rparams.b_A) == 0.05 -# assert isinstance(model.rparams.eta_g, Var) -# assert value(model.rparams.eta_g) == 0.8 -# assert isinstance(model.rparams.eta_h, Var) -# assert value(model.rparams.eta_h) == 0.8 -# assert isinstance(model.rparams.k_h, Var) -# assert value(model.rparams.k_h) == 3 -# assert isinstance(model.rparams.K_X, Var) -# assert value(model.rparams.K_X) == 0.1 -# assert isinstance(model.rparams.K_NH, Var) -# assert value(model.rparams.K_NH) == 1e-3 -# assert isinstance(model.rparams.k_a, Var) -# assert value(model.rparams.k_a) == 50 -# -# -# class TestReactionBlock(object): -# @pytest.fixture(scope="class") -# def model(self): -# model = ConcreteModel() -# model.pparams = ASM3ParameterBlock() -# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) -# -# model.props = model.pparams.build_state_block([1]) -# -# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) -# -# return model -# -# @pytest.mark.unit -# def test_build(self, model): -# assert model.rxns[1].conc_mass_comp_ref is model.props[1].conc_mass_comp -# -# @pytest.mark.unit -# def test_rxn_rate(self, model): -# assert isinstance(model.rxns[1].reaction_rate, Var) -# assert len(model.rxns[1].reaction_rate) == 8 -# assert isinstance(model.rxns[1].rate_expression, Constraint) -# assert len(model.rxns[1].rate_expression) == 8 -# -# @pytest.mark.unit -# def test_get_reaction_rate_basis(self, model): -# assert model.rxns[1].get_reaction_rate_basis() == MaterialFlowBasis.mass -# -# @pytest.mark.component -# def test_initialize(self, model): -# assert model.rxns.initialize() is None -# -# @pytest.mark.unit -# def check_units(self, model): -# assert_units_consistent(model) -# -# -# class TestASM3ReactionScaler(object): -# @pytest.mark.unit -# def test_variable_scaling_routine(self): -# model = ConcreteModel() -# model.pparams = ASM3ParameterBlock() -# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) -# -# model.props = model.pparams.build_state_block([1]) -# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) -# -# # Trigger build of reaction properties -# model.rxns[1].reaction_rate -# -# scaler = model.rxns[1].default_scaler() -# assert isinstance(scaler, ASM3ReactionScaler) -# -# scaler.variable_scaling_routine(model.rxns[1]) -# -# assert isinstance(model.rxns[1].scaling_factor, Suffix) -# -# sfx = model.rxns[1].scaling_factor -# assert len(sfx) == 8 -# assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8) -# -# @pytest.mark.unit -# def test_constraint_scaling_routine(self): -# model = ConcreteModel() -# model.pparams = ASM3ParameterBlock() -# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) -# -# model.props = model.pparams.build_state_block([1]) -# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) -# -# # Trigger build of reaction properties -# model.rxns[1].reaction_rate -# -# scaler = model.rxns[1].default_scaler() -# assert isinstance(scaler, ASM3ReactionScaler) -# -# scaler.constraint_scaling_routine(model.rxns[1]) -# -# assert isinstance(model.rxns[1].scaling_factor, Suffix) -# -# sfx = model.rxns[1].scaling_factor -# assert len(sfx) == 8 -# assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx( -# 2.380752e5, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx( -# 1.49540985e8, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx( -# 1.75226112e6, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx( -# 2.88e6, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx( -# 1.728e7, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx( -# 1.728e5, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx( -# 3.174336e5, rel=1e-8 -# ) -# assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1, rel=1e-8) -# -# @pytest.mark.unit -# def test_scale_model(self): -# model = ConcreteModel() -# model.pparams = ASM3ParameterBlock() -# model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) -# -# model.props = model.pparams.build_state_block([1]) -# model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) -# -# # Trigger build of reaction properties -# model.rxns[1].reaction_rate -# -# scaler = model.rxns[1].default_scaler() -# assert isinstance(scaler, ASM3ReactionScaler) -# -# scaler.scale_model(model.rxns[1]) -# -# assert isinstance(model.rxns[1].scaling_factor, Suffix) -# -# sfx = model.rxns[1].scaling_factor -# assert len(sfx) == 16 -# assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8) -# -# assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx(1e2, rel=1e-8) -# assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1e2, rel=1e-8) -# -# -# class TestReactor: -# @pytest.fixture(scope="class") -# def model(self): -# m = ConcreteModel() -# -# m.fs = FlowsheetBlock(dynamic=False) -# -# m.fs.props = ASM3ParameterBlock() -# m.fs.rxn_props = ASM3ReactionParameterBlock(property_package=m.fs.props) -# -# m.fs.R1 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) -# -# # Feed conditions based on manual mass balance of inlet and recycle streams -# m.fs.R1.inlet.flow_vol.fix(92230 * units.m**3 / units.day) -# m.fs.R1.inlet.temperature.fix(298.15 * units.K) -# m.fs.R1.inlet.pressure.fix(1 * units.atm) -# m.fs.R1.inlet.conc_mass_comp[0, "S_I"].fix(30 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "S_S"].fix(14.6112 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "X_I"].fix(1149 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "X_S"].fix(89.324 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "X_BH"].fix(2542.03 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "X_BA"].fix(148.6 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "X_P"].fix(448 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "S_O"].fix(0.3928 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "S_NO"].fix(8.32 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "S_NH"].fix(7.696 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "S_ND"].fix(1.9404 * units.g / units.m**3) -# m.fs.R1.inlet.conc_mass_comp[0, "X_ND"].fix(5.616 * units.g / units.m**3) -# m.fs.R1.inlet.alkalinity.fix(4.704 * units.mol / units.m**3) -# -# m.fs.R1.volume.fix(1000 * units.m**3) -# -# return m -# -# @pytest.mark.unit -# def test_dof(self, model): -# assert degrees_of_freedom(model) == 0 -# -# @pytest.mark.unit -# def test_unit_consistency(self, model): -# assert_units_consistent(model) == 0 -# -# @pytest.mark.component -# def test_solve(self, model): -# model.fs.R1.initialize() -# -# solver = get_solver() -# results = solver.solve(model, tee=True) -# assert check_optimal_termination(results) -# -# @pytest.mark.component -# def test_solution(self, model): -# assert value(model.fs.R1.outlet.flow_vol[0]) == pytest.approx(1.0675, rel=1e-4) -# assert value(model.fs.R1.outlet.temperature[0]) == pytest.approx( -# 298.15, rel=1e-4 -# ) -# assert value(model.fs.R1.outlet.pressure[0]) == pytest.approx(101325, rel=1e-4) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_I"]) == pytest.approx( -# 30e-3, rel=1e-5 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_S"]) == pytest.approx( -# 2.81e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_I"]) == pytest.approx( -# 1149e-3, rel=1e-3 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_S"]) == pytest.approx( -# 82.1e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_BH"]) == pytest.approx( -# 2552e-3, rel=1e-3 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_BA"]) == pytest.approx( -# 149e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_P"]) == pytest.approx( -# 449e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_O"]) == pytest.approx( -# 4.3e-6, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NO"]) == pytest.approx( -# 5.36e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NH"]) == pytest.approx( -# 7.92e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_ND"]) == pytest.approx( -# 1.22e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_ND"]) == pytest.approx( -# 5.29e-3, rel=1e-2 -# ) -# assert value(model.fs.R1.outlet.alkalinity[0]) == pytest.approx( -# 4.93e-3, rel=1e-2 -# ) + @pytest.mark.unit + def test_build(self, model): + assert model.rparams.reaction_block_class is ASM3ReactionBlock + + assert len(model.rparams.rate_reaction_idx) == 12 + for i in model.rparams.rate_reaction_idx: + assert i in [ + "R1", + "R2", + "R3", + "R4", + "R5", + "R6", + "R7", + "R8", + "R9", + "R10", + "R11", + "R12", + ] + + assert len(model.rparams.rate_reaction_stoichiometry) == 12 * 14 + for i in model.rparams.rate_reaction_stoichiometry: + assert i[0] in [ + "R1", + "R2", + "R3", + "R4", + "R5", + "R6", + "R7", + "R8", + "R9", + "R10", + "R11", + "R12", + ] + assert i[1] == "Liq" + assert i[2] in [ + "H2O", + "S_O", + "S_I", + "S_S", + "X_S", + "S_NH4", + "S_N2", + "S_NOX", + "S_ALK", + "X_I", + "X_S", + "X_H", + "X_STO", + "X_A", + "X_TSS", + ] + + assert isinstance(model.rparams.Y_A, Var) + assert value(model.rparams.Y_A) == 0.24 + assert isinstance(model.rparams.Y_H_O2, Var) + assert value(model.rparams.Y_H_O2) == 0.63 + assert isinstance(model.rparams.Y_H_NOX, Var) + assert value(model.rparams.Y_H_NOX) == 0.54 + assert isinstance(model.rparams.Y_STO_O2, Var) + assert value(model.rparams.Y_STO_O2) == 0.85 + assert isinstance(model.rparams.Y_STO_NOX, Var) + assert value(model.rparams.Y_STO_NOX) == 0.80 + + assert isinstance(model.rparams.k_H, Var) + assert value(model.rparams.k_H["10C"]) == 2 + assert value(model.rparams.k_H["20C"]) == 3 + assert isinstance(model.rparams.K_X, Var) + assert value(model.rparams.K_X) == 1 + assert isinstance(model.rparams.k_STO, Var) + assert value(model.rparams.k_STO["10C"]) == 2.5 + assert value(model.rparams.k_STO["20C"]) == 5 + assert isinstance(model.rparams.eta_NOX, Var) + assert value(model.rparams.eta_NOX) == 0.6 + assert isinstance(model.rparams.K_O2, Var) + assert value(model.rparams.K_O2) == 0.2e-3 + assert isinstance(model.rparams.K_NOX, Var) + assert value(model.rparams.K_NOX) == 0.5e-3 + assert isinstance(model.rparams.K_S, Var) + assert value(model.rparams.K_S) == 2e-3 + assert isinstance(model.rparams.K_STO, Var) + assert value(model.rparams.K_STO) == 1 + assert isinstance(model.rparams.mu_H, Var) + assert value(model.rparams.mu_H["10C"]) == 1 + assert value(model.rparams.mu_H["20C"]) == 2 + assert isinstance(model.rparams.K_NH4, Var) + assert value(model.rparams.K_NH4) == 0.01e-3 + assert isinstance(model.rparams.K_ALK, Var) + assert value(model.rparams.K_ALK) == 0.1e-3 + assert isinstance(model.rparams.b_H_O2, Var) + assert value(model.rparams.b_H_O2["10C"]) == 0.1 + assert value(model.rparams.b_H_O2["20C"]) == 0.2 + assert isinstance(model.rparams.b_H_NOX, Var) + assert value(model.rparams.b_H_NOX["10C"]) == 0.05 + assert value(model.rparams.b_H_NOX["20C"]) == 0.1 + assert isinstance(model.rparams.b_STO_O2, Var) + assert value(model.rparams.b_STO_O2["10C"]) == 0.1 + assert value(model.rparams.b_STO_O2["20C"]) == 0.2 + assert isinstance(model.rparams.b_STO_NOX, Var) + assert value(model.rparams.b_STO_NOX["10C"]) == 0.05 + assert value(model.rparams.b_STO_NOX["20C"]) == 0.1 + + assert isinstance(model.rparams.mu_A, Var) + assert value(model.rparams.mu_A["10C"]) == 0.35 + assert value(model.rparams.mu_A["20C"]) == 1 + assert isinstance(model.rparams.K_A_NH4, Var) + assert value(model.rparams.K_A_NH4) == 1e-3 + assert isinstance(model.rparams.K_A_O2, Var) + assert value(model.rparams.K_A_O2) == 0.5e-3 + assert isinstance(model.rparams.K_A_ALK, Var) + assert value(model.rparams.K_A_ALK) == 0.5e-3 + assert isinstance(model.rparams.b_A_O2, Var) + assert value(model.rparams.b_A_O2["10C"]) == 0.05 + assert value(model.rparams.b_A_O2["20C"]) == 0.15 + assert isinstance(model.rparams.b_A_NOX, Var) + assert value(model.rparams.b_A_NOX["10C"]) == 0.02 + assert value(model.rparams.b_A_NOX["20C"]) == 0.05 + + +class TestReactionBlock(object): + @pytest.fixture(scope="class") + def model(self): + model = ConcreteModel() + model.pparams = ASM3ParameterBlock() + model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) + + model.props = model.pparams.build_state_block([1]) + + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + return model + + @pytest.mark.unit + def test_build(self, model): + assert model.rxns[1].conc_mass_comp_ref is model.props[1].conc_mass_comp + + @pytest.mark.unit + def test_rxn_rate(self, model): + assert isinstance(model.rxns[1].reaction_rate, Var) + assert len(model.rxns[1].reaction_rate) == 12 + assert isinstance(model.rxns[1].rate_expression, Constraint) + assert len(model.rxns[1].rate_expression) == 12 + + @pytest.mark.unit + def test_get_reaction_rate_basis(self, model): + assert model.rxns[1].get_reaction_rate_basis() == MaterialFlowBasis.mass + + @pytest.mark.component + def test_initialize(self, model): + assert model.rxns.initialize() is None + + @pytest.mark.unit + def check_units(self, model): + assert_units_consistent(model) + + +class TestASM3ReactionScaler(object): + @pytest.mark.unit + def test_variable_scaling_routine(self): + model = ConcreteModel() + model.pparams = ASM3ParameterBlock() + model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) + + model.props = model.pparams.build_state_block([1]) + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + # Trigger build of reaction properties + model.rxns[1].reaction_rate + + scaler = model.rxns[1].default_scaler() + assert isinstance(scaler, ASM3ReactionScaler) + + scaler.variable_scaling_routine(model.rxns[1]) + + assert isinstance(model.rxns[1].scaling_factor, Suffix) + + sfx = model.rxns[1].scaling_factor + assert len(sfx) == 12 + assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R9"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R10"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R11"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R12"]] == pytest.approx(1e2, rel=1e-8) + + @pytest.mark.unit + def test_constraint_scaling_routine(self): + model = ConcreteModel() + model.pparams = ASM3ParameterBlock() + model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) + + model.props = model.pparams.build_state_block([1]) + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + # Trigger build of reaction properties + model.rxns[1].reaction_rate + + scaler = model.rxns[1].default_scaler() + assert isinstance(scaler, ASM3ReactionScaler) + + scaler.constraint_scaling_routine(model.rxns[1]) + + assert isinstance(model.rxns[1].scaling_factor, Suffix) + + sfx = model.rxns[1].scaling_factor + assert len(sfx) == 12 + assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx( + 5.76e5, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx( + 1.76608e5, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx( + 147909628.8, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx( + 865901.1, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx( + 721584295.2, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx( + 4328640, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx( + 4350283200, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx( + 4328640, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R9"]] == pytest.approx( + 4350283200, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R10"]] == pytest.approx( + 877441.7, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R11"]] == pytest.approx( + 5788800, rel=1e-5 + ) + assert sfx[model.rxns[1].rate_expression["R12"]] == pytest.approx( + 3490646400, rel=1e-5 + ) + + @pytest.mark.unit + def test_scale_model(self): + model = ConcreteModel() + model.pparams = ASM3ParameterBlock() + model.rparams = ASM3ReactionParameterBlock(property_package=model.pparams) + + model.props = model.pparams.build_state_block([1]) + model.rxns = model.rparams.build_reaction_block([1], state_block=model.props) + + # Trigger build of reaction properties + model.rxns[1].reaction_rate + + scaler = model.rxns[1].default_scaler() + assert isinstance(scaler, ASM3ReactionScaler) + + scaler.scale_model(model.rxns[1]) + + assert isinstance(model.rxns[1].scaling_factor, Suffix) + # + sfx = model.rxns[1].scaling_factor + assert len(sfx) == 24 + assert sfx[model.rxns[1].reaction_rate["R1"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R2"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R3"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R4"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R5"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R6"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R7"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R8"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R9"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R10"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R11"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].reaction_rate["R12"]] == pytest.approx(1e2, rel=1e-8) + + assert sfx[model.rxns[1].rate_expression["R1"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R2"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R3"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R4"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R7"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R8"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R9"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R10"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R11"]] == pytest.approx(1e2, rel=1e-8) + assert sfx[model.rxns[1].rate_expression["R12"]] == pytest.approx(1e2, rel=1e-8) + + +class TestReactor: + @pytest.fixture(scope="class") + def model(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ASM3ParameterBlock() + m.fs.rxn_props = ASM3ReactionParameterBlock(property_package=m.fs.props) + + m.fs.R1 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) + + m.fs.R1.inlet.flow_vol.fix(92230 * units.m**3 / units.day) + m.fs.R1.inlet.temperature.fix(288.15 * units.K) + m.fs.R1.inlet.pressure.fix(1 * units.atm) + m.fs.R1.inlet.conc_mass_comp[0, "S_O"].fix( + 0.0333140769653528 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "S_I"].fix(30 * units.g / units.m**3) + m.fs.R1.inlet.conc_mass_comp[0, "S_S"].fix( + 1.79253833233150 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "S_NH4"].fix( + 7.47840572528914 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "S_N2"].fix( + 25.0222401125193 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "S_NOX"].fix( + 4.49343937121928 * units.g / units.m**3 + ) + m.fs.R1.inlet.alkalinity.fix(4.95892616814772 * units.mol / units.m**3) + m.fs.R1.inlet.conc_mass_comp[0, "X_I"].fix( + 1460.88032984731 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "X_S"].fix( + 239.049918909639 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "X_H"].fix( + 1624.51533042293 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "X_STO"].fix( + 316.937373308996 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "X_A"].fix( + 130.798830163795 * units.g / units.m**3 + ) + m.fs.R1.inlet.conc_mass_comp[0, "X_TSS"].fix( + 3044.89285508125 * units.g / units.m**3 + ) + + m.fs.R1.volume.fix(1000 * units.m**3) + + return m + + @pytest.mark.unit + def test_dof(self, model): + assert degrees_of_freedom(model) == 0 + + @pytest.mark.unit + def test_unit_consistency(self, model): + assert_units_consistent(model) == 0 + + @pytest.mark.component + def test_solve(self, model): + model.fs.R1.initialize() + + solver = get_solver() + results = solver.solve(model, tee=True) + assert check_optimal_termination(results) + + @pytest.mark.component + def test_solution(self, model): + assert value(model.fs.R1.outlet.flow_vol[0]) == pytest.approx(1.0675, rel=1e-4) + assert value(model.fs.R1.outlet.temperature[0]) == pytest.approx( + 288.15, rel=1e-4 + ) + assert value(model.fs.R1.outlet.pressure[0]) == pytest.approx(101325, rel=1e-4) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_O"]) == pytest.approx( + 3.76884e-7, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_I"]) == pytest.approx( + 30e-3, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_S"]) == pytest.approx( + 4.43187e-4, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NH4"]) == pytest.approx( + 7.6437e-3, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_N2"]) == pytest.approx( + 2.7107e-2, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NOX"]) == pytest.approx( + 2.4126e-3, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_I"]) == pytest.approx( + 1.4612, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_S"]) == pytest.approx( + 0.23243, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_H"]) == pytest.approx( + 1.6264, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_STO"]) == pytest.approx( + 0.31677, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_A"]) == pytest.approx( + 0.13074, rel=1e-4 + ) + assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_TSS"]) == pytest.approx( + 3.0417, rel=1e-4 + ) + assert value(model.fs.R1.outlet.alkalinity[0]) == pytest.approx( + 4.96155e-3, rel=1e-4 + ) From 220d0c65794d5980e4e10abc04a54c7a4c8e78ed Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Mon, 16 Jun 2025 10:24:26 -0400 Subject: [PATCH 11/21] revise test --- .../unit_specific/activated_sludge/tests/test_asm3_thermo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py index 75477d8f3a..2dcf822b6b 100644 --- a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py @@ -192,7 +192,7 @@ def test_build(self, model): assert str(model.props[1].material_flow_expression[j].expr) == str( model.props[1].flow_vol * model.props[1].alkalinity - * (12 * units.kg / units.kmol) + * (61 * units.kg / units.kmol) ) else: assert str(model.props[1].material_flow_expression[j].expr) == str( From d5135cc982dd2a04edbf11b2c461166f9e93f6ed Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Mon, 16 Jun 2025 11:07:34 -0400 Subject: [PATCH 12/21] code linting --- .../unit_specific/activated_sludge/asm3_reactions.py | 1 - 1 file changed, 1 deletion(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py index 2217dde9f5..cb196a4419 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py @@ -36,7 +36,6 @@ from idaes.core.util.exceptions import BurntToast import idaes.logger as idaeslog from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme -import idaes.core.util.scaling as iscale # Some more information about this module __author__ = "Chenyu Wang, Adam Atia" From d8fa5547e5b81abc6972a4a4c69c040bd60f334c Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Mon, 16 Jun 2025 19:03:14 -0400 Subject: [PATCH 13/21] fix tests; switch cstr import to watertap cstr which includes hrt; change asr volumes; eliminate unnecessary code in ASR flowsheet --- .../activated_sludge_reactor_train.py | 84 ++++--------------- .../activated_sludge/asm1_properties.py | 4 + watertap/unit_models/cstr.py | 12 ++- .../unit_models/tests/test_aeration_tank.py | 4 +- .../unit_models/tests/test_dewatering_unit.py | 4 +- 5 files changed, 27 insertions(+), 81 deletions(-) diff --git a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py b/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py index c5be8b4032..24dfd3efec 100644 --- a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py +++ b/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py @@ -41,7 +41,6 @@ from idaes.core import FlowsheetBlock from idaes.models.unit_models import ( - CSTR, Feed, Mixer, Separator, @@ -60,7 +59,7 @@ ) from idaes.core.util.initialization import propagate_state -from watertap.unit_models import AerationTank +from watertap.unit_models import AerationTank, CSTR from watertap.property_models.unit_specific.activated_sludge.asm1_properties import ( ASM1ParameterBlock, @@ -120,26 +119,11 @@ def build_flowsheet(): m.fs.S1 = Separator( property_package=m.fs.props, outlet_list=["effluent","M1_inlet", "R2_inlet"] ) - # Clarifier - # # TODO: Replace with more detailed model when available - # m.fs.CL1 = Separator( - # property_package=m.fs.props, - # outlet_list=["underflow", "effluent"], - # split_basis=SplittingType.componentFlow, - # ) - # # Sludge purge splitter - # m.fs.SP6 = Separator( - # property_package=m.fs.props, - # outlet_list=["recycle", "waste"], - # split_basis=SplittingType.totalFlow, - # ) - # # Mixing sludge recycle and R5 underflow - # m.fs.MX6 = Mixer(property_package=m.fs.props, inlet_list=["clarifier", "reactor"]) + # Product Blocks m.fs.Treated = Product(property_package=m.fs.props) - # Recycle pressure changer - use a simple isothermal unit for now - # m.fs.P1 = PressureChanger(property_package=m.fs.props) + # Link units m.fs.feed_to_m1 = Arc(source=m.fs.feed.outlet, destination=m.fs.M1.feed) @@ -153,8 +137,7 @@ def build_flowsheet(): m.fs.s1_to_effluent = Arc(source=m.fs.S1.effluent, destination=m.fs.Treated.inlet) m.fs.s1_to_m1 = Arc(source=m.fs.S1.M1_inlet, destination=m.fs.M1.recycle) m.fs.s1_to_r2 = Arc(source=m.fs.S1.R2_inlet, destination=m.fs.R2.inlet) - # m.fs.stream14 = Arc(source=m.fs.MX6.outlet, destination=m.fs.P1.inlet) - # m.fs.stream15 = Arc(source=m.fs.P1.outlet, destination=m.fs.MX1.recycle) + pyo.TransformationFactory("network.expand_arcs").apply_to(m) return m @@ -181,9 +164,9 @@ def set_operating_conditions(m): # Reactor sizing m.fs.R1.volume.fix(1000 * pyo.units.m**3) m.fs.R2.volume.fix(1000 * pyo.units.m**3) - m.fs.R3.volume.fix(1333 * pyo.units.m**3) - m.fs.R4.volume.fix(1333 * pyo.units.m**3) - m.fs.R5.volume.fix(1333 * pyo.units.m**3) + m.fs.R3.volume.fix(1000 * pyo.units.m**3) + m.fs.R4.volume.fix(1000 * pyo.units.m**3) + m.fs.R5.volume.fix(1000 * pyo.units.m**3) # Injection rates to Reactors 2 and 4 for j in m.fs.props.component_list: @@ -191,9 +174,9 @@ def set_operating_conditions(m): # All components except S_O have no injection m.fs.R2.injection[:, :, j].fix(0) m.fs.R4.injection[:, :, j].fix(0) - # Then set injections rates for O2 - m.fs.R2.outlet.conc_mass_comp[:, "S_O"].fix(1.72e-3) - m.fs.R4.outlet.conc_mass_comp[:, "S_O"].fix(2.43e-3) + + m.fs.R2.KLa.fix(240) + m.fs.R4.KLa.fix(240) # Set fraction of outflow from reactor 5 that goes to recycle m.fs.S1.split_fraction[:, "M1_inlet"].fix(0.5) @@ -201,8 +184,6 @@ def set_operating_conditions(m): # Set fraction of outflow from reactor 5 that goes to effluent m.fs.S1.split_fraction[:, "effluent"].fix(0.4) - # # Outlet pressure from recycle pump - # m.fs.P1.outlet.pressure.fix(101325) # Check degrees of freedom print("DOF = ", degrees_of_freedom(m)) @@ -307,31 +288,7 @@ def initialize_flowsheet(m): # for o in order: # print(o[0].name) - # # Initial guesses for flow into first reactor - # tear_guesses = { - # "flow_vol": {0: 1.0675}, - # "conc_mass_comp": { - # (0, "S_I"): 0.03, - # (0, "S_S"): 0.0146, - # (0, "X_I"): 1.149, - # (0, "X_S"): 0.0893, - # (0, "X_BH"): 2.542, - # (0, "X_BA"): 0.149, - # (0, "X_P"): 0.448, - # (0, "S_O"): 3.93e-4, - # (0, "S_NO"): 8.32e-3, - # (0, "S_NH"): 7.7e-3, - # (0, "S_ND"): 1.94e-3, - # (0, "X_ND"): 5.62e-3, - # }, - # "alkalinity": {0: 5e-3}, - # "temperature": {0: 298.15}, - # "pressure": {0: 101325}, - # } - - # # Pass the tear_guess to the SD tool - # seq.set_guesses_for(m.fs.R1.inlet, tear_guesses) - + def function(unit): unit.initialize(outlvl=idaeslog.DEBUG) @@ -340,23 +297,8 @@ def function(unit): def solve_flowsheet(m): # Solve overall flowsheet to close recycle loop solver = get_solver() - results = solver.solve(m) - check_solve(results, checkpoint="closing recycle", logger=_log, fail_flag=True) - - # Switch to fixed KLa in R3 and R4 (S_O concentration is controlled in R5) - m.fs.R2.KLa.fix(10) - m.fs.R4.KLa.fix(10) - m.fs.R2.outlet.conc_mass_comp[:, "S_O"].unfix() - m.fs.R4.outlet.conc_mass_comp[:, "S_O"].unfix() - - # Resolve with controls in place results = solver.solve(m, tee=True) - check_solve( - results, - checkpoint="re-solve with controls in place", - logger=_log, - fail_flag=True, - ) + check_solve(results, checkpoint="closing recycle", logger=_log, fail_flag=True) return m, results @@ -394,3 +336,5 @@ def solve_flowsheet(m): time_point=0, ) print(stream_table_dataframe_to_string(stream_table)) + m.fs.R2._get_performance_contents() + m.fs.R4._get_performance_contents() \ No newline at end of file diff --git a/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py index d91dff088f..4117332f29 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm1_properties.py @@ -533,6 +533,10 @@ def define_state_vars(self): def define_display_vars(self): return { + "TSS": self.TSS, + "COD": self.COD, + "BOD5": self.BOD5, + "TKN": self.TKN, "Volumetric Flowrate": self.flow_vol, "Molar Alkalinity": self.alkalinity, "Mass Concentration": self.conc_mass_comp, diff --git a/watertap/unit_models/cstr.py b/watertap/unit_models/cstr.py index d77ebb8b7d..8e66c40fa6 100644 --- a/watertap/unit_models/cstr.py +++ b/watertap/unit_models/cstr.py @@ -66,18 +66,16 @@ def build(self): doc="Hydraulic retention time", ) - def CSTR_retention_time_rule(self, t): + @self.Constraint( + self.flowsheet().time, + doc="Hydraulic retention time equation", + ) + def eq_hydraulic_retention_time(self, t): return ( self.hydraulic_retention_time[t] == self.volume[t] / self.control_volume.properties_in[t].flow_vol ) - self.CSTR_retention_time = Constraint( - self.flowsheet().time, - rule=CSTR_retention_time_rule, - doc="Total CSTR retention time", - ) - @property def default_costing_method(self): return cost_cstr diff --git a/watertap/unit_models/tests/test_aeration_tank.py b/watertap/unit_models/tests/test_aeration_tank.py index 343a5d961e..52a3fb9e63 100644 --- a/watertap/unit_models/tests/test_aeration_tank.py +++ b/watertap/unit_models/tests/test_aeration_tank.py @@ -175,9 +175,9 @@ def test_build(self, model): == model.fs.unit.control_volume.mass_transfer_term[k] ) - assert number_variables(model) == 100 + assert number_variables(model) == 99 assert number_total_constraints(model) == 49 - assert number_unused_variables(model) == 1 + assert number_unused_variables(model) == 0 @pytest.mark.component def test_units(self, model): diff --git a/watertap/unit_models/tests/test_dewatering_unit.py b/watertap/unit_models/tests/test_dewatering_unit.py index 50a157ae62..e2a7244f5d 100644 --- a/watertap/unit_models/tests/test_dewatering_unit.py +++ b/watertap/unit_models/tests/test_dewatering_unit.py @@ -188,9 +188,9 @@ def test_build(self, du): assert hasattr(du.fs.unit.overflow, "pressure") assert hasattr(du.fs.unit.overflow, "alkalinity") - assert number_variables(du) == 83 + assert number_variables(du) == 82 assert number_total_constraints(du) == 62 - assert number_unused_variables(du) == 4 + assert number_unused_variables(du) == 3 @pytest.mark.unit def test_dof(self, du): From 5aa6f63043baf3707b5cc8ebbef58af80dd4b39a Mon Sep 17 00:00:00 2001 From: Chenyu Date: Wed, 25 Jun 2025 13:41:39 -0400 Subject: [PATCH 14/21] add metrics --- .../activated_sludge/asm3_properties.py | 160 +++++++++--------- 1 file changed, 76 insertions(+), 84 deletions(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py index c1867d9fee..0b5a3e28b8 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py @@ -223,13 +223,13 @@ def build(self): # domain=pyo.PositiveReals, # doc="Conversion factor applied for TSS calculation", # ) - # self.BOD5_factor = pyo.Param( - # ["raw", "effluent"], - # initialize={"raw": 0.65, "effluent": 0.25}, - # units=pyo.units.dimensionless, - # domain=pyo.PositiveReals, - # doc="Conversion factor for BOD5", - # ) + self.BOD5_factor = pyo.Param( + ["raw", "effluent"], + initialize={"raw": 0.65, "effluent": 0.25}, + units=pyo.units.dimensionless, + domain=pyo.PositiveReals, + doc="Conversion factor for BOD5", + ) # Fix Vars that are treated as Params for v in self.component_objects(pyo.Var): v.fix() @@ -247,11 +247,11 @@ def define_metadata(cls, obj): obj.define_custom_properties( { "alkalinity": {"method": None}, - # "TSS": {"method": "_TSS"}, + "TSS": {"method": "_TSS"}, # "BOD5": {"method": "_BOD5"}, - # "TKN": {"method": "_TKN"}, - # "Total_N": {"method": "_Total_N"}, - # "COD": {"method": "_COD"}, + "TKN": {"method": "_TKN"}, + "Total_N": {"method": "_Total_N"}, + "COD": {"method": "_COD"}, } ) obj.add_default_units( @@ -503,79 +503,71 @@ def energy_density_expression(self): rule=energy_density_expression, doc="Energy density term" ) - # def _TSS(self): - # tss = ( - # self.conc_mass_comp["X_S"] - # + self.conc_mass_comp["X_I"] - # + self.conc_mass_comp["X_BH"] - # + self.conc_mass_comp["X_BA"] - # + self.conc_mass_comp["X_P"] - # ) - # return self.params.COD_to_SS * tss - # - # self.TSS = pyo.Expression( - # rule=_TSS, - # doc="Total suspended solids (TSS)", - # ) - # - # def _BOD5(self, i): - # bod5 = ( - # self.conc_mass_comp["S_S"] - # + self.conc_mass_comp["X_S"] - # + (1 - self.params.f_p) - # * (self.conc_mass_comp["X_BH"] + self.conc_mass_comp["X_BA"]) - # ) - # # TODO: 0.25 should be a parameter instead as it changes by influent/effluent - # return self.params.BOD5_factor[i] * bod5 - # - # self.BOD5 = pyo.Expression( - # ["raw", "effluent"], - # rule=_BOD5, - # doc="Five-day Biological Oxygen Demand (BOD5)", - # ) - # - # def _COD(self): - # cod = ( - # self.conc_mass_comp["S_S"] - # + self.conc_mass_comp["S_I"] - # + self.conc_mass_comp["X_S"] - # + self.conc_mass_comp["X_I"] - # + self.conc_mass_comp["X_BH"] - # + self.conc_mass_comp["X_BA"] - # + self.conc_mass_comp["X_P"] - # ) - # return cod - # - # self.COD = pyo.Expression( - # rule=_COD, - # doc="Chemical Oxygen Demand", - # ) - # - # def _TKN(self): - # tkn = ( - # self.conc_mass_comp["S_NH4"] - # + self.conc_mass_comp["S_ND"] - # + self.conc_mass_comp["X_ND"] - # + self.params.i_xb - # * (self.conc_mass_comp["X_BH"] + self.conc_mass_comp["X_BA"]) - # + self.params.i_xp - # * (self.conc_mass_comp["X_P"] + self.conc_mass_comp["X_I"]) - # ) - # return tkn - # - # self.TKN = pyo.Expression( - # rule=_TKN, - # doc="Total Kjeldahl Nitrogen", - # ) - # - # def _Total_N(self): - # totaln = self.TKN + self.conc_mass_comp["S_NOX"] - # return totaln - # - # self.Total_N = pyo.Expression( - # rule=_Total_N, - # doc="Total Nitrogen", - # ) + def _TSS(self): + tss = self.conc_mass_comp["X_TSS"] + return tss + + self.TSS = pyo.Expression( + rule=_TSS, + doc="Total suspended solids (TSS)", + ) + + def _BOD5(self, i): + bod5 = ( + self.conc_mass_comp["S_S"] + + (1 - self.params.f_SI) * self.conc_mass_comp["X_S"] + + (1 - self.params.f_XIH) * self.conc_mass_comp["X_H"] + + (1 - self.params.f_XIA) * self.conc_mass_comp["X_AUT"] + ) + # TODO: 0.25 should be a parameter instead as it changes by influent/effluent + return self.params.BOD5_factor[i] * bod5 + + self.BOD5 = pyo.Expression( + ["raw", "effluent"], + rule=_BOD5, + doc="Five-day Biological Oxygen Demand (BOD5)", + ) + + def _COD(self): + cod = ( + self.conc_mass_comp["S_S"] + + self.conc_mass_comp["S_I"] + + self.conc_mass_comp["X_S"] + + self.conc_mass_comp["X_I"] + + self.conc_mass_comp["X_STO"] + ) + return cod + + self.COD = pyo.Expression( + rule=_COD, + doc="Chemical Oxygen Demand", + ) + + def _TKN(self): + tkn = ( + self.conc_mass_comp["S_NH4"] + + self.params.i_NSS * self.conc_mass_comp["S_S"] + + self.params.i_NSI * self.conc_mass_comp["S_I"] + + self.params.i_NXI * self.conc_mass_comp["X_I"] + + self.params.i_NXS * self.conc_mass_comp["X_S"] + + self.params.i_NBM + * (self.conc_mass_comp["X_H"] + self.conc_mass_comp["X_A"]) + ) + return tkn + + self.TKN = pyo.Expression( + rule=_TKN, + doc="Total Kjeldahl Nitrogen", + ) + + def _Total_N(self): + totaln = self.TKN + self.conc_mass_comp["S_NOX"] + return totaln + + self.Total_N = pyo.Expression( + rule=_Total_N, + doc="Total Nitrogen", + ) def get_material_flow_terms(self, p, j): return self.material_flow_expression[j] From 957f997c7c925e5a019664415d55a394558e2ad2 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Wed, 25 Jun 2025 16:07:45 -0400 Subject: [PATCH 15/21] delete BOD5 --- .../activated_sludge/asm3_properties.py | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py index 0b5a3e28b8..fcd302c0b4 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py @@ -223,13 +223,13 @@ def build(self): # domain=pyo.PositiveReals, # doc="Conversion factor applied for TSS calculation", # ) - self.BOD5_factor = pyo.Param( - ["raw", "effluent"], - initialize={"raw": 0.65, "effluent": 0.25}, - units=pyo.units.dimensionless, - domain=pyo.PositiveReals, - doc="Conversion factor for BOD5", - ) + # self.BOD5_factor = pyo.Param( + # ["raw", "effluent"], + # initialize={"raw": 0.65, "effluent": 0.25}, + # units=pyo.units.dimensionless, + # domain=pyo.PositiveReals, + # doc="Conversion factor for BOD5", + # ) # Fix Vars that are treated as Params for v in self.component_objects(pyo.Var): v.fix() @@ -512,21 +512,21 @@ def _TSS(self): doc="Total suspended solids (TSS)", ) - def _BOD5(self, i): - bod5 = ( - self.conc_mass_comp["S_S"] - + (1 - self.params.f_SI) * self.conc_mass_comp["X_S"] - + (1 - self.params.f_XIH) * self.conc_mass_comp["X_H"] - + (1 - self.params.f_XIA) * self.conc_mass_comp["X_AUT"] - ) - # TODO: 0.25 should be a parameter instead as it changes by influent/effluent - return self.params.BOD5_factor[i] * bod5 - - self.BOD5 = pyo.Expression( - ["raw", "effluent"], - rule=_BOD5, - doc="Five-day Biological Oxygen Demand (BOD5)", - ) + # def _BOD5(self, i): + # bod5 = ( + # self.conc_mass_comp["S_S"] + # + (1 - self.params.f_SI) * self.conc_mass_comp["X_S"] + # + (1 - self.params.f_XIH) * self.conc_mass_comp["X_H"] + # + (1 - self.params.f_XIA) * self.conc_mass_comp["X_AUT"] + # ) + # # TODO: 0.25 should be a parameter instead as it changes by influent/effluent + # return self.params.BOD5_factor[i] * bod5 + # + # self.BOD5 = pyo.Expression( + # ["raw", "effluent"], + # rule=_BOD5, + # doc="Five-day Biological Oxygen Demand (BOD5)", + # ) def _COD(self): cod = ( From d9f8d76d4ae5387cd8a9727cef03758037496cb5 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Wed, 25 Jun 2025 16:11:28 -0400 Subject: [PATCH 16/21] add metrics and test --- .../activated_sludge/asm3_properties.py | 48 +------------------ .../tests/test_asm3_thermo.py | 15 +++--- 2 files changed, 7 insertions(+), 56 deletions(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py index fcd302c0b4..6a8b3e6e8b 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_properties.py @@ -199,37 +199,7 @@ def build(self): domain=pyo.PositiveReals, doc="TSS to COD ratio for X_STO based on PHB (g-TSS / g-X_STO)", ) - # self.f_p = pyo.Var( - # initialize=0.08, - # units=pyo.units.dimensionless, - # domain=pyo.PositiveReals, - # doc="Fraction of biomass yielding particulate products, f_p", - # ) - # self.i_xb = pyo.Var( - # initialize=0.08, - # units=pyo.units.dimensionless, - # domain=pyo.PositiveReals, - # doc="Mass fraction of N per COD in biomass, i_xb", - # ) - # self.i_xp = pyo.Var( - # initialize=0.06, - # units=pyo.units.dimensionless, - # domain=pyo.PositiveReals, - # doc="Mass fraction of N per COD in particulates, i_xp", - # ) - # self.COD_to_SS = pyo.Var( - # initialize=0.75, - # units=pyo.units.dimensionless, - # domain=pyo.PositiveReals, - # doc="Conversion factor applied for TSS calculation", - # ) - # self.BOD5_factor = pyo.Param( - # ["raw", "effluent"], - # initialize={"raw": 0.65, "effluent": 0.25}, - # units=pyo.units.dimensionless, - # domain=pyo.PositiveReals, - # doc="Conversion factor for BOD5", - # ) + # Fix Vars that are treated as Params for v in self.component_objects(pyo.Var): v.fix() @@ -512,22 +482,6 @@ def _TSS(self): doc="Total suspended solids (TSS)", ) - # def _BOD5(self, i): - # bod5 = ( - # self.conc_mass_comp["S_S"] - # + (1 - self.params.f_SI) * self.conc_mass_comp["X_S"] - # + (1 - self.params.f_XIH) * self.conc_mass_comp["X_H"] - # + (1 - self.params.f_XIA) * self.conc_mass_comp["X_AUT"] - # ) - # # TODO: 0.25 should be a parameter instead as it changes by influent/effluent - # return self.params.BOD5_factor[i] * bod5 - # - # self.BOD5 = pyo.Expression( - # ["raw", "effluent"], - # rule=_BOD5, - # doc="Five-day Biological Oxygen Demand (BOD5)", - # ) - def _COD(self): cod = ( self.conc_mass_comp["S_S"] diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py index 2dcf822b6b..b360ad7eda 100644 --- a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_thermo.py @@ -341,15 +341,12 @@ def test_initialize(self, model): def check_units(self, model): assert_units_consistent(model) - -# @pytest.mark.unit -# def test_expressions(self, model): -# assert value(model.props[1].TSS) == 0.375 -# assert value(model.props[1].COD) == pytest.approx(0.7, rel=1e-3) -# assert value(model.props[1].BOD5["effluent"]) == 0.096 -# assert value(model.props[1].BOD5["raw"]) == 0.096 * 0.65 / 0.25 -# assert value(model.props[1].TKN) == pytest.approx(0.328, rel=1e-3) -# assert value(model.props[1].Total_N) == pytest.approx(0.428, rel=1e-3) + @pytest.mark.unit + def test_expressions(self, model): + assert value(model.props[1].TSS) == 0.1 + assert value(model.props[1].COD) == pytest.approx(0.5, rel=1e-3) + assert value(model.props[1].TKN) == pytest.approx(0.124, rel=1e-3) + assert value(model.props[1].Total_N) == pytest.approx(0.224, rel=1e-3) class TestASM3PropertiesScaler: From 11c1762e944fbbfc28e3025c3c244d5d6135783c Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Thu, 26 Jun 2025 15:44:30 -0400 Subject: [PATCH 17/21] rename file after adding asm3 compatibility --- ..._sludge_reactor_train.py => UConn_WRRF.py} | 167 +++++++++++------- 1 file changed, 106 insertions(+), 61 deletions(-) rename watertap/flowsheets/activated_sludge/{activated_sludge_reactor_train.py => UConn_WRRF.py} (71%) diff --git a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py b/watertap/flowsheets/activated_sludge/UConn_WRRF.py similarity index 71% rename from watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py rename to watertap/flowsheets/activated_sludge/UConn_WRRF.py index 24dfd3efec..c3c2cd6182 100644 --- a/watertap/flowsheets/activated_sludge/activated_sludge_reactor_train.py +++ b/watertap/flowsheets/activated_sludge/UConn_WRRF.py @@ -10,8 +10,9 @@ # "https://github.com/watertap-org/watertap/" ################################################################################# """ -Example of activated sludge process model. -Initial layout: +Example of activated sludge process model, based on 5 CSTR representation of the University of Connecticut's Water Resource Recovery Facility, conceptualized by Stuber's Process Systems and Operations Research Laboratory. + +Layout: * 5 reactors * R1, R3, R5: anoxic * R2 and R4: aerobic @@ -27,27 +28,16 @@ * Anoxic reactors as standard CSTRs (CSTR) * Aerobic reactors as AerationTanks -Note that a pressure changer is required in the recycle stream to ensure the -pressure inside the recycle loop is bounded. As the inlet Mixer uses a pressure -minimization constraint and there is no pressure drop in the reactors, if pressure -is not specified at some point within the recycle loop then it becomes unbounded. - """ __author__ = "Adam Atia" +from enum import Enum, auto import pyomo.environ as pyo from pyomo.network import Arc, SequentialDecomposition from idaes.core import FlowsheetBlock -from idaes.models.unit_models import ( - Feed, - Mixer, - Separator, - PressureChanger, - Product, - MomentumMixingType -) +from idaes.models.unit_models import Feed, Mixer, Separator, Product, MomentumMixingType from idaes.models.unit_models.separator import SplittingType from watertap.core.solvers import get_solver from idaes.core.util.model_statistics import degrees_of_freedom @@ -67,64 +57,72 @@ from watertap.property_models.unit_specific.activated_sludge.asm1_reactions import ( ASM1ReactionParameterBlock, ) - +from watertap.property_models.unit_specific.activated_sludge.asm3_properties import ( + ASM3ParameterBlock, +) +from watertap.property_models.unit_specific.activated_sludge.asm3_reactions import ( + ASM3ReactionParameterBlock, +) from watertap.core.util.initialization import check_solve, interval_initializer # Set up logger _log = idaeslog.getLogger(__name__) -def build_flowsheet(): +class ASMModel(auto): + asm1 = auto() + asm3 = auto() + + +def build_flowsheet(asm_model=ASMModel.asm1): m = pyo.ConcreteModel() m.fs = FlowsheetBlock(dynamic=False) + if asm_model == ASMModel.asm1: + m.fs.props = ASM1ParameterBlock() + m.fs.rxn_props = ASM1ReactionParameterBlock(property_package=m.fs.props) + elif asm_model == ASMModel.asm3: + m.fs.props = ASM3ParameterBlock() + m.fs.rxn_props = ASM3ReactionParameterBlock(property_package=m.fs.props) - m.fs.props = ASM1ParameterBlock() - m.fs.rxn_props = ASM1ReactionParameterBlock(property_package=m.fs.props) # Feed water stream m.fs.feed = Feed(property_package=m.fs.props) # Mixer for feed water and recycled sludge - m.fs.M1 = Mixer(property_package=m.fs.props, - inlet_list=["feed", "recycle"], - momentum_mixing_type=MomentumMixingType.none - ) + m.fs.M1 = Mixer( + property_package=m.fs.props, + inlet_list=["feed", "recycle"], + momentum_mixing_type=MomentumMixingType.none, + ) m.fs.M1.outlet.pressure.fix() # m.fs.M1.pressure_equality_constraints[0,2].deactivate() - # First reactor (anoxic) + # First reactor (anoxic) m.fs.R1 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) - - # Second reactor (aerobic) + + # Second reactor (aerobic) m.fs.R2 = AerationTank(property_package=m.fs.props, reaction_package=m.fs.rxn_props) - - + # Mixer for R1 and R2 outlets - m.fs.M2= Mixer(property_package=m.fs.props, inlet_list=["R1_outlet", "R2_outlet"], - momentum_mixing_type=MomentumMixingType.none - ) + m.fs.M2 = Mixer( + property_package=m.fs.props, + inlet_list=["R1_outlet", "R2_outlet"], + momentum_mixing_type=MomentumMixingType.none, + ) m.fs.M2.outlet.pressure.fix() # Third reactor (anoxic) - m.fs.R3 = CSTR( - property_package=m.fs.props, reaction_package=m.fs.rxn_props - ) + m.fs.R3 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) # Fourth reactor (aerobic) - CSTR with injection - m.fs.R4 = AerationTank( - property_package=m.fs.props, reaction_package=m.fs.rxn_props - ) - # Fifth reactor (anoxic) - m.fs.R5 = CSTR( - property_package=m.fs.props, reaction_package=m.fs.rxn_props - ) + m.fs.R4 = AerationTank(property_package=m.fs.props, reaction_package=m.fs.rxn_props) + # Fifth reactor (anoxic) + m.fs.R5 = CSTR(property_package=m.fs.props, reaction_package=m.fs.rxn_props) m.fs.S1 = Separator( - property_package=m.fs.props, outlet_list=["effluent","M1_inlet", "R2_inlet"] + property_package=m.fs.props, outlet_list=["effluent", "M1_inlet", "R2_inlet"] ) - # Product Blocks m.fs.Treated = Product(property_package=m.fs.props) - # Link units m.fs.feed_to_m1 = Arc(source=m.fs.feed.outlet, destination=m.fs.M1.feed) m.fs.m1_to_r1 = Arc(source=m.fs.M1.outlet, destination=m.fs.R1.inlet) @@ -139,11 +137,11 @@ def build_flowsheet(): m.fs.s1_to_r2 = Arc(source=m.fs.S1.R2_inlet, destination=m.fs.R2.inlet) pyo.TransformationFactory("network.expand_arcs").apply_to(m) - + return m -def set_operating_conditions(m): - # Feed Water Conditions + +def set_asm1_inlet_conditions(m): m.fs.feed.flow_vol.fix(18446 * pyo.units.m**3 / pyo.units.day) m.fs.feed.temperature.fix(298.15 * pyo.units.K) m.fs.feed.pressure.fix(1 * pyo.units.atm) @@ -161,6 +159,55 @@ def set_operating_conditions(m): m.fs.feed.conc_mass_comp[0, "X_ND"].fix(10.59 * pyo.units.g / pyo.units.m**3) m.fs.feed.alkalinity.fix(7 * pyo.units.mol / pyo.units.m**3) + +def set_asm3_inlet_conditions(m): + m.fs.feed.flow_vol.fix(92230 * pyo.units.m**3 / pyo.units.day) + m.fs.feed.temperature.fix(288.15 * pyo.units.K) + m.fs.feed.pressure.fix(1 * pyo.units.atm) + m.fs.feed.conc_mass_comp[0, "S_O"].fix( + 0.0333140769653528 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "S_I"].fix(30 * pyo.units.g / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "S_S"].fix( + 1.79253833233150 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "S_NH4"].fix( + 7.47840572528914 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "S_N2"].fix( + 25.0222401125193 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "S_NOX"].fix( + 4.49343937121928 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.alkalinity.fix(4.95892616814772 * pyo.units.mol / pyo.units.m**3) + m.fs.feed.conc_mass_comp[0, "X_I"].fix( + 1460.88032984731 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "X_S"].fix( + 239.049918909639 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "X_H"].fix( + 1624.51533042293 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "X_STO"].fix( + 316.937373308996 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "X_A"].fix( + 130.798830163795 * pyo.units.g / pyo.units.m**3 + ) + m.fs.feed.conc_mass_comp[0, "X_TSS"].fix( + 3044.89285508125 * pyo.units.g / pyo.units.m**3 + ) + + +def set_operating_conditions(m, asm_model=ASMModel.asm1): + # Feed Water Conditions + if asm_model == ASMModel.asm1: + set_asm1_inlet_conditions(m) + elif asm_model == ASMModel.asm3: + set_asm3_inlet_conditions(m) + # Reactor sizing m.fs.R1.volume.fix(1000 * pyo.units.m**3) m.fs.R2.volume.fix(1000 * pyo.units.m**3) @@ -184,11 +231,11 @@ def set_operating_conditions(m): # Set fraction of outflow from reactor 5 that goes to effluent m.fs.S1.split_fraction[:, "effluent"].fix(0.4) - # Check degrees of freedom print("DOF = ", degrees_of_freedom(m)) assert degrees_of_freedom(m) == 0 + def scale_flowsheet(m): # Apply scaling for var in m.fs.component_data_objects(pyo.Var, descend_into=True): @@ -202,12 +249,13 @@ def scale_flowsheet(m): iscale.set_scaling_factor(var, 1e2) iscale.calculate_scaling_factors(m.fs) -def init_and_propagate(blk,arc=None,source=None,destination=None): + +def init_and_propagate(blk, arc=None, source=None, destination=None): blk.initialize() if arc is not None: propagate_state(arc) elif source is not None: - propagate_state(source=source,destination=destination) + propagate_state(source=source, destination=destination) def initialize_flowsheet(m): @@ -227,7 +275,7 @@ def initialize_flowsheet(m): m.fs.M2.initialize() propagate_state(m.fs.m2_to_r3) - m.fs.R3.initialize() + m.fs.R3.initialize() propagate_state(m.fs.r3_to_r4) m.fs.R4.initialize() @@ -253,7 +301,7 @@ def initialize_flowsheet(m): m.fs.M2.initialize() propagate_state(m.fs.m2_to_r3) - m.fs.R3.initialize() + m.fs.R3.initialize() propagate_state(m.fs.r3_to_r4) m.fs.R4.initialize() @@ -266,12 +314,11 @@ def initialize_flowsheet(m): propagate_state(m.fs.s1_to_effluent) propagate_state(m.fs.s1_to_m1) propagate_state(m.fs.s1_to_r2) - + m.fs.Treated.initialize() interval_initializer(m) - # Apply sequential decomposition - 1 iteration should suffice seq = SequentialDecomposition() seq.options.select_tear_method = "heuristic" @@ -288,12 +335,12 @@ def initialize_flowsheet(m): # for o in order: # print(o[0].name) - def function(unit): unit.initialize(outlvl=idaeslog.DEBUG) seq.run(m, function) + def solve_flowsheet(m): # Solve overall flowsheet to close recycle loop solver = get_solver() @@ -307,8 +354,8 @@ def solve_flowsheet(m): # This method builds and runs a steady state activated sludge # flowsheet. # m, results = build_flowsheet() - m = build_flowsheet() - set_operating_conditions(m) + m = build_flowsheet(asm_model=ASMModel.asm3) + set_operating_conditions(m, asm_model=ASMModel.asm3) scale_flowsheet(m) initialize_flowsheet(m) @@ -317,7 +364,6 @@ def solve_flowsheet(m): m, res = solve_flowsheet(m) - stream_table = create_stream_table_dataframe( { "Feed": m.fs.feed.outlet, @@ -330,11 +376,10 @@ def solve_flowsheet(m): "R5": m.fs.R5.outlet, "S1 to M1": m.fs.S1.M1_inlet, "S1 to R2": m.fs.S1.R2_inlet, - "Effluent": m.fs.Treated.inlet - + "Effluent": m.fs.Treated.inlet, }, time_point=0, ) print(stream_table_dataframe_to_string(stream_table)) m.fs.R2._get_performance_contents() - m.fs.R4._get_performance_contents() \ No newline at end of file + m.fs.R4._get_performance_contents() From 6ff8190ebe86b4e142f6cb3357759a353f9d0464 Mon Sep 17 00:00:00 2001 From: luohezhiming Date: Mon, 30 Jun 2025 16:24:35 -0400 Subject: [PATCH 18/21] revise kinetic equation for R5 --- .../activated_sludge/asm3_reactions.py | 4 ++++ .../tests/test_asm3_reaction.py | 20 +++++++++---------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py index cb196a4419..be93386a69 100644 --- a/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py +++ b/watertap/property_models/unit_specific/activated_sludge/asm3_reactions.py @@ -723,6 +723,10 @@ def rate_expression_rule(b, r): b.params.K_O2 / (b.params.K_O2 + b.conc_mass_comp_ref["S_O"]) ) + * ( + b.conc_mass_comp_ref["S_NOX"] + / (b.params.K_NOX + b.conc_mass_comp_ref["S_NOX"]) + ) * ( b.conc_mass_comp_ref["S_NH4"] / (b.params.K_NH4 + b.conc_mass_comp_ref["S_NH4"]) diff --git a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py index bf1850e08c..6e0fa91195 100644 --- a/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py +++ b/watertap/property_models/unit_specific/activated_sludge/tests/test_asm3_reaction.py @@ -288,7 +288,7 @@ def test_constraint_scaling_routine(self): 865901.1, rel=1e-5 ) assert sfx[model.rxns[1].rate_expression["R5"]] == pytest.approx( - 721584295.2, rel=1e-5 + 725192216.7, rel=1e-5 ) assert sfx[model.rxns[1].rate_expression["R6"]] == pytest.approx( 4328640, rel=1e-5 @@ -439,22 +439,22 @@ def test_solution(self, model): ) assert value(model.fs.R1.outlet.pressure[0]) == pytest.approx(101325, rel=1e-4) assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_O"]) == pytest.approx( - 3.76884e-7, rel=1e-4 + 3.77072e-7, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_I"]) == pytest.approx( 30e-3, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_S"]) == pytest.approx( - 4.43187e-4, rel=1e-4 + 4.38103e-4, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NH4"]) == pytest.approx( - 7.6437e-3, rel=1e-4 + 7.6828e-3, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_N2"]) == pytest.approx( - 2.7107e-2, rel=1e-4 + 2.69511e-2, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "S_NOX"]) == pytest.approx( - 2.4126e-3, rel=1e-4 + 2.56815e-3, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_I"]) == pytest.approx( 1.4612, rel=1e-4 @@ -463,17 +463,17 @@ def test_solution(self, model): 0.23243, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_H"]) == pytest.approx( - 1.6264, rel=1e-4 + 1.6259, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_STO"]) == pytest.approx( - 0.31677, rel=1e-4 + 0.31777, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_A"]) == pytest.approx( 0.13074, rel=1e-4 ) assert value(model.fs.R1.outlet.conc_mass_comp[0, "X_TSS"]) == pytest.approx( - 3.0417, rel=1e-4 + 3.04183, rel=1e-4 ) assert value(model.fs.R1.outlet.alkalinity[0]) == pytest.approx( - 4.96155e-3, rel=1e-4 + 4.9614e-3, rel=1e-4 ) From 6b95b8dee9cb7282bb65508527c073160e91e7ee Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Mon, 30 Jun 2025 16:44:24 -0400 Subject: [PATCH 19/21] start test file --- .../flowsheets/activated_sludge/UConn_WRRF.py | 4 +- .../activated_sludge/tests/test_UConn_WRRF.py | 143 ++++++++++++++++++ 2 files changed, 145 insertions(+), 2 deletions(-) create mode 100644 watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py diff --git a/watertap/flowsheets/activated_sludge/UConn_WRRF.py b/watertap/flowsheets/activated_sludge/UConn_WRRF.py index c3c2cd6182..2edb01e16d 100644 --- a/watertap/flowsheets/activated_sludge/UConn_WRRF.py +++ b/watertap/flowsheets/activated_sludge/UConn_WRRF.py @@ -347,7 +347,7 @@ def solve_flowsheet(m): results = solver.solve(m, tee=True) check_solve(results, checkpoint="closing recycle", logger=_log, fail_flag=True) - return m, results + return results if __name__ == "__main__": @@ -362,7 +362,7 @@ def solve_flowsheet(m): scale_flowsheet(m) - m, res = solve_flowsheet(m) + res = solve_flowsheet(m) stream_table = create_stream_table_dataframe( { diff --git a/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py b/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py new file mode 100644 index 0000000000..6cf754ed56 --- /dev/null +++ b/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py @@ -0,0 +1,143 @@ +################################################################################# +# 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/" +################################################################################# +""" +Tests for UConn WRRF flowsheet example. + +""" + +__author__ = "Adam Atia" + +import pytest + +from pyomo.environ import assert_optimal_termination, value +from pyomo.util.check_units import assert_units_consistent + +from idaes.core.util.model_statistics import degrees_of_freedom + +from watertap.flowsheets.activated_sludge import UConn_WRRF as uconn +from watertap.flowsheets.activated_sludge.UConn_WRRF import ASMModel + + + +class TestUConnFlowsheet: + @pytest.fixture(scope="class") + def model(self): + m = uconn.build_flowsheet(asm_model=ASMModel.asm3) + + uconn.set_operating_conditions(m, asm_model=ASMModel.asm3) + uconn.scale_flowsheet(m) + + uconn.initialize_flowsheet(m) + + uconn.scale_flowsheet(m) + + res = uconn.solve_flowsheet(m) + + m.results = res + + return m + + @pytest.mark.integration + def test_structure(self, model): + assert_units_consistent(model) + assert degrees_of_freedom(model) == 0 + assert_optimal_termination(model.results) + + @pytest.mark.integration + def test_results(self, model): + # Treated water + assert value(model.fs.Treated.flow_vol[0]) == pytest.approx(0.20904, rel=1e-4) + assert value(model.fs.Treated.temperature[0]) == pytest.approx(298.15, rel=1e-4) + assert value(model.fs.Treated.pressure[0]) == pytest.approx(101325, rel=1e-4) + assert value(model.fs.Treated.conc_mass_comp[0, "S_I"]) == pytest.approx( + 30e-3, rel=1e-5 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "S_S"]) == pytest.approx( + 8.89e-4, rel=1e-2 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "X_I"]) == pytest.approx( + 4.39e-3, rel=1e-3 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "X_S"]) == pytest.approx( + 1.88e-4, rel=1e-2 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "X_BH"]) == pytest.approx( + 9.78e-3, rel=1e-3 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "X_BA"]) == pytest.approx( + 5.73e-4, rel=1e-2 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "X_P"]) == pytest.approx( + 1.73e-3, rel=1e-2 + ) + # S_O is slightly off, but probably due to differences in injection rates + assert value(model.fs.Treated.conc_mass_comp[0, "S_O"]) == pytest.approx( + 4.49e-4, rel=1e-2 + ) + # Slightly off in last significant digit + assert value(model.fs.Treated.conc_mass_comp[0, "S_NO"]) == pytest.approx( + 10.14e-3, rel=1e-2 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "S_NH"]) == pytest.approx( + 1.86e-3, rel=1e-2 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "S_ND"]) == pytest.approx( + 6.88e-4, rel=1e-2 + ) + assert value(model.fs.Treated.conc_mass_comp[0, "X_ND"]) == pytest.approx( + 1.35e-5, rel=1e-2 + ) + assert value(model.fs.Treated.alkalinity[0]) == pytest.approx(4.13e-3, rel=1e-2) + + # Sludge stream + assert value(model.fs.Sludge.flow_vol[0]) == pytest.approx(4.457e-3, rel=1e-4) + assert value(model.fs.Sludge.temperature[0]) == pytest.approx(298.15, rel=1e-4) + assert value(model.fs.Sludge.pressure[0]) == pytest.approx(101325, rel=1e-4) + assert value(model.fs.Sludge.conc_mass_comp[0, "S_I"]) == pytest.approx( + 30e-3, rel=1e-5 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "S_S"]) == pytest.approx( + 8.89e-4, rel=1e-2 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "X_I"]) == pytest.approx( + 2247e-3, rel=1e-3 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "X_S"]) == pytest.approx( + 96.8e-3, rel=1e-2 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "X_BH"]) == pytest.approx( + 5004e-3, rel=1e-3 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "X_BA"]) == pytest.approx( + 292e-3, rel=1e-2 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "X_P"]) == pytest.approx( + 884e-3, rel=1e-2 + ) + # S_O is slightly off, but probably due to differences in injection rates + assert value(model.fs.Sludge.conc_mass_comp[0, "S_O"]) == pytest.approx( + 4.49e-4, rel=1e-2 + ) + # Slightly off in last significant digit + assert value(model.fs.Sludge.conc_mass_comp[0, "S_NO"]) == pytest.approx( + 10.14e-3, rel=1e-2 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "S_NH"]) == pytest.approx( + 1.86e-3, rel=1e-2 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "S_ND"]) == pytest.approx( + 6.88e-4, rel=1e-2 + ) + assert value(model.fs.Sludge.conc_mass_comp[0, "X_ND"]) == pytest.approx( + 6.92e-3, rel=1e-2 + ) + assert value(model.fs.Sludge.alkalinity[0]) == pytest.approx(4.13e-3, rel=1e-2) From f5bd9e98ce24c14e6d0b7bfb900fda68bbf796e9 Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Tue, 1 Jul 2025 16:09:49 -0400 Subject: [PATCH 20/21] add some testing, update seq decomp to direct instead of wegstein, run black --- .../flowsheets/activated_sludge/UConn_WRRF.py | 2 +- .../activated_sludge/tests/test_UConn_WRRF.py | 132 +++++++----------- 2 files changed, 53 insertions(+), 81 deletions(-) diff --git a/watertap/flowsheets/activated_sludge/UConn_WRRF.py b/watertap/flowsheets/activated_sludge/UConn_WRRF.py index 2edb01e16d..ea59afa5e6 100644 --- a/watertap/flowsheets/activated_sludge/UConn_WRRF.py +++ b/watertap/flowsheets/activated_sludge/UConn_WRRF.py @@ -322,7 +322,7 @@ def initialize_flowsheet(m): # Apply sequential decomposition - 1 iteration should suffice seq = SequentialDecomposition() seq.options.select_tear_method = "heuristic" - seq.options.tear_method = "Wegstein" + seq.options.tear_method = "Direct" seq.options.iterLim = 1 G = seq.create_graph(m) diff --git a/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py b/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py index 6cf754ed56..8cf74e659a 100644 --- a/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py +++ b/watertap/flowsheets/activated_sludge/tests/test_UConn_WRRF.py @@ -23,26 +23,31 @@ from idaes.core.util.model_statistics import degrees_of_freedom -from watertap.flowsheets.activated_sludge import UConn_WRRF as uconn +from watertap.flowsheets.activated_sludge.UConn_WRRF import ( + build_flowsheet, + set_operating_conditions, + scale_flowsheet, + initialize_flowsheet, + solve_flowsheet, +) from watertap.flowsheets.activated_sludge.UConn_WRRF import ASMModel - -class TestUConnFlowsheet: +class TestUConnFlowsheetASM3: @pytest.fixture(scope="class") def model(self): - m = uconn.build_flowsheet(asm_model=ASMModel.asm3) + m = build_flowsheet(asm_model=ASMModel.asm3) - uconn.set_operating_conditions(m, asm_model=ASMModel.asm3) - uconn.scale_flowsheet(m) + set_operating_conditions(m, asm_model=ASMModel.asm3) + scale_flowsheet(m) - uconn.initialize_flowsheet(m) + initialize_flowsheet(m) - uconn.scale_flowsheet(m) + scale_flowsheet(m) - res = uconn.solve_flowsheet(m) + # res = solve_flowsheet(m) - m.results = res + # m.results = res return m @@ -50,94 +55,61 @@ def model(self): def test_structure(self, model): assert_units_consistent(model) assert degrees_of_freedom(model) == 0 - assert_optimal_termination(model.results) + + @pytest.mark.integration + def test_solve(self, model): + res = solve_flowsheet(model) + assert_optimal_termination(res) @pytest.mark.integration def test_results(self, model): + comps = model.fs.props.solute_set + model.fs.Treated.display() + # Treated water - assert value(model.fs.Treated.flow_vol[0]) == pytest.approx(0.20904, rel=1e-4) - assert value(model.fs.Treated.temperature[0]) == pytest.approx(298.15, rel=1e-4) + assert value(model.fs.Treated.flow_vol[0]) == pytest.approx( + 1.06747685185185, rel=1e-4 + ) + assert value(model.fs.Treated.temperature[0]) == pytest.approx(288.15, rel=1e-4) assert value(model.fs.Treated.pressure[0]) == pytest.approx(101325, rel=1e-4) assert value(model.fs.Treated.conc_mass_comp[0, "S_I"]) == pytest.approx( 30e-3, rel=1e-5 ) - assert value(model.fs.Treated.conc_mass_comp[0, "S_S"]) == pytest.approx( - 8.89e-4, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "S_N2"]) == pytest.approx( + 0.02888, rel=1e-2 ) - assert value(model.fs.Treated.conc_mass_comp[0, "X_I"]) == pytest.approx( - 4.39e-3, rel=1e-3 + assert value(model.fs.Treated.conc_mass_comp[0, "S_NH4"]) == pytest.approx( + 0.001519, rel=1e-2 ) - assert value(model.fs.Treated.conc_mass_comp[0, "X_S"]) == pytest.approx( - 1.88e-4, rel=1e-2 - ) - assert value(model.fs.Treated.conc_mass_comp[0, "X_BH"]) == pytest.approx( - 9.78e-3, rel=1e-3 - ) - assert value(model.fs.Treated.conc_mass_comp[0, "X_BA"]) == pytest.approx( - 5.73e-4, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "S_NOX"]) == pytest.approx( + 0.007174, rel=1e-2 ) - assert value(model.fs.Treated.conc_mass_comp[0, "X_P"]) == pytest.approx( - 1.73e-3, rel=1e-2 - ) - # S_O is slightly off, but probably due to differences in injection rates assert value(model.fs.Treated.conc_mass_comp[0, "S_O"]) == pytest.approx( - 4.49e-4, rel=1e-2 - ) - # Slightly off in last significant digit - assert value(model.fs.Treated.conc_mass_comp[0, "S_NO"]) == pytest.approx( - 10.14e-3, rel=1e-2 - ) - assert value(model.fs.Treated.conc_mass_comp[0, "S_NH"]) == pytest.approx( - 1.86e-3, rel=1e-2 + 0.00099, rel=1e-2 ) - assert value(model.fs.Treated.conc_mass_comp[0, "S_ND"]) == pytest.approx( - 6.88e-4, rel=1e-2 - ) - assert value(model.fs.Treated.conc_mass_comp[0, "X_ND"]) == pytest.approx( - 1.35e-5, rel=1e-2 - ) - assert value(model.fs.Treated.alkalinity[0]) == pytest.approx(4.13e-3, rel=1e-2) - - # Sludge stream - assert value(model.fs.Sludge.flow_vol[0]) == pytest.approx(4.457e-3, rel=1e-4) - assert value(model.fs.Sludge.temperature[0]) == pytest.approx(298.15, rel=1e-4) - assert value(model.fs.Sludge.pressure[0]) == pytest.approx(101325, rel=1e-4) - assert value(model.fs.Sludge.conc_mass_comp[0, "S_I"]) == pytest.approx( - 30e-3, rel=1e-5 - ) - assert value(model.fs.Sludge.conc_mass_comp[0, "S_S"]) == pytest.approx( - 8.89e-4, rel=1e-2 - ) - assert value(model.fs.Sludge.conc_mass_comp[0, "X_I"]) == pytest.approx( - 2247e-3, rel=1e-3 - ) - assert value(model.fs.Sludge.conc_mass_comp[0, "X_S"]) == pytest.approx( - 96.8e-3, rel=1e-2 - ) - assert value(model.fs.Sludge.conc_mass_comp[0, "X_BH"]) == pytest.approx( - 5004e-3, rel=1e-3 + assert value(model.fs.Treated.conc_mass_comp[0, "S_S"]) == pytest.approx( + 0.00016, rel=1e-2 ) - assert value(model.fs.Sludge.conc_mass_comp[0, "X_BA"]) == pytest.approx( - 292e-3, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_A"]) == pytest.approx( + 0.13165, rel=1e-2 ) - assert value(model.fs.Sludge.conc_mass_comp[0, "X_P"]) == pytest.approx( - 884e-3, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_H"]) == pytest.approx( + 1.632358, rel=1e-2 ) - # S_O is slightly off, but probably due to differences in injection rates - assert value(model.fs.Sludge.conc_mass_comp[0, "S_O"]) == pytest.approx( - 4.49e-4, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_I"]) == pytest.approx( + 1.46378, rel=1e-2 ) - # Slightly off in last significant digit - assert value(model.fs.Sludge.conc_mass_comp[0, "S_NO"]) == pytest.approx( - 10.14e-3, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_S"]) == pytest.approx( + 0.209026, rel=1e-2 ) - assert value(model.fs.Sludge.conc_mass_comp[0, "S_NH"]) == pytest.approx( - 1.86e-3, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_S"]) == pytest.approx( + 0.209026, rel=1e-2 ) - assert value(model.fs.Sludge.conc_mass_comp[0, "S_ND"]) == pytest.approx( - 6.88e-4, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_STO"]) == pytest.approx( + 0.304689, rel=1e-2 ) - assert value(model.fs.Sludge.conc_mass_comp[0, "X_ND"]) == pytest.approx( - 6.92e-3, rel=1e-2 + assert value(model.fs.Treated.conc_mass_comp[0, "X_TSS"]) == pytest.approx( + 3.02503, rel=1e-3 ) - assert value(model.fs.Sludge.alkalinity[0]) == pytest.approx(4.13e-3, rel=1e-2) + + assert value(model.fs.Treated.alkalinity[0]) == pytest.approx(0.00495, rel=1e-2) From 52b13e36446cceec1382ecf67339e405dea46afc Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Wed, 2 Jul 2025 20:51:37 -0400 Subject: [PATCH 21/21] attempting to validate [wip] --- .../flowsheets/activated_sludge/UConn_WRRF.py | 75 ++++++++++++++++++- 1 file changed, 71 insertions(+), 4 deletions(-) diff --git a/watertap/flowsheets/activated_sludge/UConn_WRRF.py b/watertap/flowsheets/activated_sludge/UConn_WRRF.py index ea59afa5e6..2dd70c9022 100644 --- a/watertap/flowsheets/activated_sludge/UConn_WRRF.py +++ b/watertap/flowsheets/activated_sludge/UConn_WRRF.py @@ -222,11 +222,11 @@ def set_operating_conditions(m, asm_model=ASMModel.asm1): m.fs.R2.injection[:, :, j].fix(0) m.fs.R4.injection[:, :, j].fix(0) - m.fs.R2.KLa.fix(240) - m.fs.R4.KLa.fix(240) + m.fs.R2.KLa.fix(10 / pyo.units.hour) + m.fs.R4.KLa.fix(10 / pyo.units.hour) - # Set fraction of outflow from reactor 5 that goes to recycle - m.fs.S1.split_fraction[:, "M1_inlet"].fix(0.5) + # Set fraction of outflow from reactor 5 that recycles to M1 mixer + m.fs.S1.split_fraction[:, "M1_inlet"].fix(0.1) # Set fraction of outflow from reactor 5 that goes to effluent m.fs.S1.split_fraction[:, "effluent"].fix(0.4) @@ -350,12 +350,49 @@ def solve_flowsheet(m): return results +def reset_asm3_inlet_conditions(m, ini_dict): + + for k in ini_dict.keys(): + if k in m.fs.props.solute_set: + m.fs.feed.conc_mass_comp[0, k].fix( + ini_dict[k] * pyo.units.g / pyo.units.m**3 + ) + elif k == "alkalinity": + m.fs.feed.alkalinity[0].fix(ini_dict[k] * pyo.units.mol / pyo.units.m**3) + elif k == "temperature": + m.fs.feed.temperature[0].fix(ini_dict[k] + 273.15) + + elif k == "flow_vol": + m.fs.feed.flow_vol[0].fix(ini_dict[k] * pyo.units.m**3 / pyo.units.day) + else: + raise + + if __name__ == "__main__": # This method builds and runs a steady state activated sludge # flowsheet. # m, results = build_flowsheet() m = build_flowsheet(asm_model=ASMModel.asm3) set_operating_conditions(m, asm_model=ASMModel.asm3) + ini2 = { + "S_O": 2.00000074088136, + "S_I": 30, + "S_S": 1.99999995166373, + "S_NH4": 20.0000000011385, + "S_N2": 1.30283567480963e-18, + "S_NOX": 9.07971541001396e-10, + "alkalinity": 5.00000000001647, + "X_I": 100.000000001382, + "X_S": 39.9999999624405, + "X_H": 100.000000012367, + "X_STO": 40.0000000397251, + "X_A": 1.0000000001811, + "X_TSS": 200.000000007995, + "flow_vol": 36892, + "temperature": 14.8581001531874, + } + + reset_asm3_inlet_conditions(m, ini2) scale_flowsheet(m) initialize_flowsheet(m) @@ -383,3 +420,33 @@ def solve_flowsheet(m): print(stream_table_dataframe_to_string(stream_table)) m.fs.R2._get_performance_contents() m.fs.R4._get_performance_contents() + + solution = { + "S_O": 6.07975, + "S_I": 30.0, + "S_S": 0.428786, + "S_NH4": 19.8235, + "S_N2": 0.0350126, + "S_NOX": 0.283642, + "alkalinity": 4.96736, + "X_I": 100.643, + "X_S": 31.7643, + "X_H": 103.536, + "X_STO": 39.5003, + "X_A": 1.08686, + "X_TSS": 197.27, + } + + my_solution = { + k: pyo.value(m.fs.Treated.conc_mass_comp[0, k]) * 1e3 + for k in solution.keys() + if k != "alkalinity" + } + import pandas as pd + + df = pd.DataFrame({"watertap": my_solution, "julia": solution}) + df.loc["alkalinity", "watertap"] = pyo.value(m.fs.Treated.alkalinity[0]) * 1e3 + df["percent_difference"] = (df["watertap"] - df["julia"]) / df["julia"] * 100 + df["percent_difference"] = df["percent_difference"].round(1) + + df