From 66eb651ff8ec256a1d0b86072753c8af2bb3a031 Mon Sep 17 00:00:00 2001 From: Adam Atia Date: Fri, 1 Dec 2023 11:15:03 -0500 Subject: [PATCH] Add Thickener Costing (#1216) * add initial costing for thickener * add testing for thickener mods and costing * modfiy bsm2 --- tutorials/BSM2.ipynb | 20 ++++ watertap/costing/unit_models/thickener.py | 64 ++++++++++++ .../BSM2.py | 4 + .../unit_models/tests/test_dewatering_unit.py | 1 - .../unit_models/tests/test_thickener_unit.py | 98 +++++++++++++++++-- watertap/unit_models/thickener.py | 92 ++++++++++++++++- 6 files changed, 267 insertions(+), 12 deletions(-) create mode 100644 watertap/costing/unit_models/thickener.py diff --git a/tutorials/BSM2.ipynb b/tutorials/BSM2.ipynb index d0c595a12c..66ad55fac8 100644 --- a/tutorials/BSM2.ipynb +++ b/tutorials/BSM2.ipynb @@ -634,6 +634,26 @@ "m.fs.DU.hydraulic_retention_time.fix(1800 * pyo.units.s)" ] }, + { + "cell_type": "markdown", + "id": "f9d6d962", + "metadata": {}, + "source": [ + "Similarly, the thickener unit includes the same equation, as well as an equation relating the thickener's dimensions. Here, we fix hydraulic retention time and thickener diameter to satisfy 0 degrees of freedom." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dfb8d9a3", + "metadata": {}, + "outputs": [], + "source": [ + "# Thickener unit\n", + "m.fs.TU.hydraulic_retention_time.fix(86400 * pyo.units.s)\n", + "m.fs.TU.diameter.fix(10 * pyo.units.m)" + ] + }, { "cell_type": "markdown", "id": "663b9e78-c1b3-41b5-b9d1-fe1fc5acd274", diff --git a/watertap/costing/unit_models/thickener.py b/watertap/costing/unit_models/thickener.py new file mode 100644 index 0000000000..a6a1ba0141 --- /dev/null +++ b/watertap/costing/unit_models/thickener.py @@ -0,0 +1,64 @@ +################################################################################# +# 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 pyomo.environ as pyo +from ..util import ( + register_costing_parameter_block, + make_capital_cost_var, +) + +""" +Ref: W. McGivney, S. Kawamura, Cost estimating manual for water treatment facilities, John Wiley & Sons, 2008. http://onlinelibrary.wiley.com/book/10.1002/9780470260036. +""" + + +def build_cost_param_block(blk): + # NOTE: costing data are for gravity sludge thickener for McGivney & Kawamura, 2008 + blk.capital_a_parameter = pyo.Var( + initialize=4729.8, + doc="A parameter for capital cost", + units=pyo.units.USD_2007 / (pyo.units.feet), + ) + blk.capital_b_parameter = pyo.Var( + initialize=37068, + doc="B parameter for capital cost", + units=pyo.units.USD_2007, + ) + + +@register_costing_parameter_block( + build_rule=build_cost_param_block, + parameter_block_name="thickener", +) +def cost_thickener(blk, cost_electricity_flow=True): + """ + Gravity Sludge Thickener costing method + """ + make_capital_cost_var(blk) + cost_blk = blk.costing_package.thickener + t0 = blk.flowsheet().time.first() + x = diameter = pyo.units.convert(blk.unit_model.diameter, to_units=pyo.units.feet) + blk.capital_cost_constraint = pyo.Constraint( + expr=blk.capital_cost + == pyo.units.convert( + cost_blk.capital_a_parameter * x + cost_blk.capital_b_parameter, + to_units=blk.costing_package.base_currency, + ) + ) + if cost_electricity_flow: + blk.costing_package.cost_flow( + pyo.units.convert( + blk.unit_model.electricity_consumption[t0], + to_units=pyo.units.kW, + ), + "electricity", + ) diff --git a/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/BSM2.py b/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/BSM2.py index 6435fc351c..ae6b0f6cc4 100644 --- a/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/BSM2.py +++ b/watertap/examples/flowsheets/case_studies/full_water_resource_recovery_facility/BSM2.py @@ -374,6 +374,10 @@ def set_operating_conditions(m): # Dewatering Unit - fix either HRT or volume. m.fs.DU.hydraulic_retention_time.fix(1800 * pyo.units.s) + # Thickener unit + m.fs.TU.hydraulic_retention_time.fix(86400 * pyo.units.s) + m.fs.TU.diameter.fix(10 * pyo.units.m) + def initialize_system(m): # Initialize flowsheet diff --git a/watertap/unit_models/tests/test_dewatering_unit.py b/watertap/unit_models/tests/test_dewatering_unit.py index 45f3f14a01..eaf8745d81 100644 --- a/watertap/unit_models/tests/test_dewatering_unit.py +++ b/watertap/unit_models/tests/test_dewatering_unit.py @@ -68,7 +68,6 @@ DewateringType, ) from idaes.core import UnitModelCostingBlock -from watertap.costing import WaterTAPCosting __author__ = "Alejandro Garciadiego, Adam Atia" diff --git a/watertap/unit_models/tests/test_thickener_unit.py b/watertap/unit_models/tests/test_thickener_unit.py index d1a2b580d9..d0074fb0e7 100644 --- a/watertap/unit_models/tests/test_thickener_unit.py +++ b/watertap/unit_models/tests/test_thickener_unit.py @@ -59,6 +59,8 @@ ModifiedASM2dParameterBlock, ) from pyomo.util.check_units import assert_units_consistent +from watertap.costing import WaterTAPCosting +from idaes.core import UnitModelCostingBlock __author__ = "Alejandro Garciadiego, Adam Atia" @@ -141,6 +143,9 @@ def tu(self): m.fs.unit.inlet.conc_mass_comp[0, "X_ND"].fix(4.7411 * units.mg / units.liter) m.fs.unit.inlet.alkalinity.fix(4.5646 * units.mol / units.m**3) + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + return m @pytest.mark.build @@ -171,8 +176,8 @@ def test_build(self, tu): assert hasattr(tu.fs.unit.overflow, "pressure") assert hasattr(tu.fs.unit.overflow, "alkalinity") - assert number_variables(tu) == 76 - assert number_total_constraints(tu) == 60 + assert number_variables(tu) == 81 + assert number_total_constraints(tu) == 63 assert number_unused_variables(tu) == 0 @pytest.mark.unit @@ -251,6 +256,9 @@ def test_solution(self, tu): assert pytest.approx(0.004564, rel=1e-3) == value( tu.fs.unit.overflow.alkalinity[0] ) + assert pytest.approx(3.82, rel=1e-3) == value(tu.fs.unit.height) + assert pytest.approx(300, rel=1e-3) == value(tu.fs.unit.volume[0]) + assert pytest.approx(78.54, rel=1e-3) == value(tu.fs.unit.surface_area) @pytest.mark.solver @pytest.mark.skipif(solver is None, reason="Solver not available") @@ -331,6 +339,9 @@ def tu_asm2d(self): ) m.fs.unit.inlet.alkalinity[0].fix(4.6663 * units.mmol / units.liter) + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + return m @pytest.mark.build @@ -361,8 +372,8 @@ def test_build(self, tu_asm2d): assert hasattr(tu_asm2d.fs.unit.overflow, "pressure") assert hasattr(tu_asm2d.fs.unit.overflow, "alkalinity") - assert number_variables(tu_asm2d) == 106 - assert number_total_constraints(tu_asm2d) == 84 + assert number_variables(tu_asm2d) == 111 + assert number_total_constraints(tu_asm2d) == 87 assert number_unused_variables(tu_asm2d) == 0 @pytest.mark.unit @@ -539,6 +550,9 @@ def tu_mod_asm2d(self): 118.3582 * units.mg / units.liter ) + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + return m @pytest.mark.build @@ -566,8 +580,8 @@ def test_build(self, tu_mod_asm2d): assert hasattr(tu_mod_asm2d.fs.unit.overflow, "temperature") assert hasattr(tu_mod_asm2d.fs.unit.overflow, "pressure") - assert number_variables(tu_mod_asm2d) == 110 - assert number_total_constraints(tu_mod_asm2d) == 83 + assert number_variables(tu_mod_asm2d) == 115 + assert number_total_constraints(tu_mod_asm2d) == 86 assert number_unused_variables(tu_mod_asm2d) == 0 @pytest.mark.unit @@ -698,3 +712,75 @@ def test_conservation(self, tu_mod_asm2d): @pytest.mark.unit def test_report(self, tu_mod_asm2d): tu_mod_asm2d.fs.unit.report() + + +@pytest.mark.solver +@pytest.mark.skipif(solver is None, reason="Solver not available") +@pytest.mark.component +def test_costing(): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = ModifiedASM2dParameterBlock() + + m.fs.unit = Thickener( + property_package=m.fs.props, + activated_sludge_model=ActivatedSludgeModelType.modified_ASM2D, + ) + + # NOTE: Concentrations of exactly 0 result in singularities, use EPS instead + EPS = 1e-8 + + m.fs.unit.inlet.flow_vol.fix(300 * units.m**3 / units.day) + m.fs.unit.inlet.temperature.fix(308.15 * units.K) + m.fs.unit.inlet.pressure.fix(1 * units.atm) + m.fs.unit.inlet.conc_mass_comp[0, "S_O2"].fix(7.9707 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_N2"].fix(29.0603 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NH4"].fix(8.0209 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_NO3"].fix(6.6395 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_PO4"].fix(7.8953 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_F"].fix(0.4748 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_A"].fix(0.0336 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_I"].fix(30 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_K"].fix(7 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_Mg"].fix(6 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "S_IC"].fix(10 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_I"].fix(1695.7695 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_S"].fix(68.2975 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_H"].fix(1855.5067 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_PAO"].fix(214.5319 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_PP"].fix(63.5316 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_PHA"].fix(2.7381 * units.mg / units.liter) + m.fs.unit.inlet.conc_mass_comp[0, "X_AUT"].fix(118.3582 * units.mg / units.liter) + + m.fs.unit.hydraulic_retention_time.fix() + m.fs.unit.diameter.fix() + + m.fs.costing = WaterTAPCosting() + + m.fs.unit.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing) + + m.fs.costing.cost_process() + + assert degrees_of_freedom(m) == 0 + + results = solver.solve(m) + + assert_optimal_termination(results) + assert_units_consistent(m) + assert hasattr(m.fs.costing, "thickener") + assert value(m.fs.costing.thickener.capital_a_parameter) == 4729.8 + assert value(m.fs.costing.thickener.capital_b_parameter) == 37068 + + # Check solutions + assert pytest.approx(220675.79, rel=1e-5) == value(m.fs.unit.costing.capital_cost) + assert pytest.approx(220675.79, rel=1e-5) == value( + units.convert( + (4729.8 * value(units.convert(10 * units.m, to_units=units.feet)) + 37068) + * units.USD_2007, + to_units=m.fs.costing.base_currency, + ) + ) + assert pytest.approx(12.5 * 0.01255, rel=1e-5) == value( + m.fs.unit.electricity_consumption[0] + ) diff --git a/watertap/unit_models/thickener.py b/watertap/unit_models/thickener.py index 4ff038103e..9d36701d74 100644 --- a/watertap/unit_models/thickener.py +++ b/watertap/unit_models/thickener.py @@ -30,12 +30,14 @@ declare_process_block_class, ) from idaes.models.unit_models.separator import SeparatorData, SplittingType - +from idaes.core.util.constants import Constants from idaes.core.util.tables import create_stream_table_dataframe import idaes.logger as idaeslog from pyomo.environ import ( Param, + Var, + NonNegativeReals, units as pyunits, ) from pyomo.common.config import ConfigValue, In @@ -43,6 +45,7 @@ from idaes.core.util.exceptions import ( ConfigurationError, ) +from watertap.costing.unit_models.thickener import cost_thickener __author__ = "Alejandro Garciadiego, Adam Atia" @@ -128,6 +131,72 @@ def build(self): doc="Fraction of suspended solids removed", ) + self.electricity_consumption = Var( + self.flowsheet().time, + units=pyunits.kW, + bounds=(0, None), + doc="Electricity consumption of unit", + ) + + # 0.01255 kWh/m3 average for centrifuge thickening in relation to flow capacity + self.energy_electric_flow_vol_inlet = Param( + self.flowsheet().time, + units=pyunits.kWh / (pyunits.m**3), + initialize=0.01255, + mutable=True, + doc="Specific electricity intensity of unit", + ) + + @self.Constraint(self.flowsheet().time, doc="Electricity consumption equation") + def eq_electricity_consumption(blk, t): + return blk.electricity_consumption[t] == pyunits.convert( + blk.energy_electric_flow_vol_inlet[t] * blk.inlet.flow_vol[t], + to_units=pyunits.kW, + ) + + self.hydraulic_retention_time = Var( + self.flowsheet().time, + initialize=86400, + domain=NonNegativeReals, + units=pyunits.s, + doc="Hydraulic retention time", + ) + self.volume = Var( + self.flowsheet().time, + initialize=1800, + domain=NonNegativeReals, + units=pyunits.m**3, + doc="Hydraulic retention time", + ) + + self.diameter = Var( + initialize=10, + domain=NonNegativeReals, + units=pyunits.m, + doc="Thickener diameter", + ) + self.height = Var( + initialize=5, + domain=NonNegativeReals, + units=pyunits.m, + doc="Thickener height", + ) + + @self.Constraint(self.flowsheet().time, doc="Hydraulic retention time equation") + def eq_hydraulic_retention(blk, t): + return ( + self.hydraulic_retention_time[t] * self.inlet.flow_vol[t] + == self.volume[t] + ) + + @self.Expression(doc="Surface area of circular thickener") + def surface_area(blk): + return self.diameter**2 * Constants.pi / 4 + + @self.Constraint(self.flowsheet().time, doc="Total volume equation") + def eq_volume(blk, t): + return self.surface_area * self.height == self.volume[t] + @self.Expression(self.flowsheet().time, doc="Suspended solids concentration") def TSS_in(blk, t): if blk.config.activated_sludge_model == ActivatedSludgeModelType.ASM1: @@ -176,14 +245,23 @@ def non_particulate_components(blk, t, i): return blk.split_fraction[t, "overflow", i] == 1 - blk.f_q_du[t] def _get_performance_contents(self, time_point=0): + var_dict = {} + expr_dict = {} + param_dict = {} if hasattr(self, "split_fraction"): - var_dict = {} for k in self.split_fraction.keys(): if k[0] == time_point: var_dict[f"Split Fraction [{str(k[1:])}]"] = self.split_fraction[k] - return {"vars": var_dict} - else: - return None + var_dict["Electricity consumption"] = self.electricity_consumption[time_point] + param_dict[ + "Specific electricity consumption" + ] = self.energy_electric_flow_vol_inlet[time_point] + var_dict["Unit Volume"] = self.volume[time_point] + var_dict["Hydraulic Retention Time"] = self.hydraulic_retention_time[time_point] + var_dict["Unit Height"] = self.height + var_dict["Unit Diameter"] = self.diameter + expr_dict["Surface Area"] = self.surface_area + return {"vars": var_dict, "params": param_dict, "exprs": expr_dict} def _get_stream_table_contents(self, time_point=0): outlet_list = self.create_outlet_list() @@ -198,3 +276,7 @@ def _get_stream_table_contents(self, time_point=0): io_dict[o] = getattr(self, o + "_state") return create_stream_table_dataframe(io_dict, time_point=time_point) + + @property + def default_costing_method(self): + return cost_thickener