Skip to content

Commit

Permalink
Add sweep analysis for electroNP flowsheet (#1077)
Browse files Browse the repository at this point in the history 
* add

* delete redundant files

* add sensitivity analysis

* add sweep analysis

* add sweep analysis

* add more analysis

* add analysis

* revise file

* add more analysis

* revise flowsheet

* revise analysis

* add sweep analysis

* delete scaling factors

* revise test file

* revise opex frac

* plinting code

* Update watertap/costing/units/electroNP.py

Co-authored-by: Adam Atia <aatia@keylogic.com>

* Update watertap/costing/units/electroNP.py

Co-authored-by: Adam Atia <aatia@keylogic.com>

* Update watertap/costing/units/electroNP.py

Co-authored-by: Adam Atia <aatia@keylogic.com>

* Update watertap/examples/flowsheets/case_studies/electroNP/multi_sweep.py

Co-authored-by: Adam Atia <aatia@keylogic.com>

* clean up the code

* revise value

* add test for display_costing

---------

Co-authored-by: Adam Atia <aatia@keylogic.com>
  • Loading branch information
@luohezhiming @adam-a-a
luohezhiming and adam-a-a authored Jul 19, 2023
1 parent d9408fd commit 06538b4
Show file tree
Hide file tree
Showing 8 changed files with 297 additions and 124 deletions.
33 changes: 31 additions & 2 deletions watertap/costing/units/electroNP.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,36 @@ def build_electroNP_cost_param_block(blk):
units=pyo.units.hr,
)
blk.sizing_cost = pyo.Var(
initialize=1.25,
initialize=1000,
doc="Reactor sizing cost",
units=pyo.units.USD_2020 / pyo.units.m**3,
)

costing = blk.parent_block()
blk.magnesium_chloride_cost = pyo.Param(
mutable=True,
initialize=0.0786,
doc="Magnesium chloride cost",
units=pyo.units.USD_2020 / pyo.units.kg,
)
costing.add_defined_flow("magnesium chloride", blk.magnesium_chloride_cost)

blk.phosphorus_recovery_value = pyo.Param(
mutable=True,
initialize=-0.07,
doc="Phosphorus recovery value",
units=pyo.units.USD_2020 / pyo.units.kg,
)
costing.add_defined_flow("phosphorus salt product", blk.phosphorus_recovery_value)


@register_costing_parameter_block(
build_rule=build_electroNP_cost_param_block,
parameter_block_name="electroNP",
)
def cost_electroNP(blk, cost_electricity_flow=True, cost_MgCl2_flow=True):
def cost_electroNP(
blk, cost_electricity_flow=True, cost_MgCl2_flow=True, cost_phosphorus_flow=True
):
"""
ElectroNP costing method
"""
Expand Down Expand Up @@ -64,6 +83,16 @@ def cost_electroNP(blk, cost_electricity_flow=True, cost_MgCl2_flow=True):
"magnesium chloride",
)

if cost_phosphorus_flow:
blk.costing_package.cost_flow(
pyo.units.convert(
blk.unit_model.byproduct.flow_vol[t0]
* blk.unit_model.byproduct.conc_mass_comp[t0, "S_PO4"],
to_units=pyo.units.kg / pyo.units.hr,
),
"phosphorus salt product",
)


def cost_electroNP_capital(blk, HRT, sizing_cost):
"""
Expand Down
8 changes: 0 additions & 8 deletions watertap/costing/watertap_costing_package.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,14 +176,6 @@ def build_global_params(self):
units=pyo.units.kg / pyo.units.kWh,
)

self.magnesium_chloride_cost = pyo.Param(
mutable=True,
initialize=0.0786,
doc="Magnesium chloride cost",
units=pyo.units.USD_2020 / pyo.units.kg,
)
self.add_defined_flow("magnesium chloride", self.magnesium_chloride_cost)

# fix the parameters
self.fix_all_vars()

Expand Down
174 changes: 77 additions & 97 deletions watertap/examples/flowsheets/case_studies/electroNP/electroNP_flowsheet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
units,
value,
assert_optimal_termination,
units as pyunits,
)
from pyomo.network import Arc
from idaes.core import (
Expand All @@ -35,18 +36,19 @@
from watertap.property_models.anaerobic_digestion.modified_adm1_reactions import (
ModifiedADM1ReactionParameterBlock,
)
from watertap.property_models.activated_sludge.simple_modified_asm2d_properties import (
SimpleModifiedASM2dParameterBlock,
from watertap.property_models.activated_sludge.modified_asm2d_properties import (
ModifiedASM2dParameterBlock,
)
from watertap.unit_models.translators.translator_adm1_simple_asm2d import (
Translator_ADM1_Simple_ASM2D,
from watertap.unit_models.translators.translator_adm1_asm2d import (
Translator_ADM1_ASM2D,
)
from watertap.unit_models.electroNP_ZO import ElectroNPZO
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.core.util.initialization import check_solve
from watertap.costing import WaterTAPCosting

# Set up logger
Expand All @@ -56,19 +58,17 @@
def build_flowsheet():
# flowsheet set up
m = pyo.ConcreteModel()

m.fs = FlowsheetBlock(dynamic=False)

m.fs.props_ADM1 = ModifiedADM1ParameterBlock()
m.fs.props_vap_ADM1 = ADM1_vaporParameterBlock()
m.fs.rxn_props_ADM1 = ModifiedADM1ReactionParameterBlock(
property_package=m.fs.props_ADM1
)
m.fs.props_ASM2D = SimpleModifiedASM2dParameterBlock(
additional_solute_list=["S_K", "S_Mg"]
)
m.fs.props_ASM2D = ModifiedASM2dParameterBlock()
m.fs.costing = WaterTAPCosting()

# Unit models
m.fs.AD = AD(
liquid_property_package=m.fs.props_ADM1,
vapor_property_package=m.fs.props_vap_ADM1,
Expand All @@ -78,7 +78,7 @@ def build_flowsheet():
)
m.fs.AD.costing = UnitModelCostingBlock(flowsheet_costing_block=m.fs.costing)

m.fs.translator_adm1_asm2d = Translator_ADM1_Simple_ASM2D(
m.fs.translator_adm1_asm2d = Translator_ADM1_ASM2D(
inlet_property_package=m.fs.props_ADM1,
outlet_property_package=m.fs.props_ASM2D,
reaction_package=m.fs.rxn_props_ADM1,
Expand All @@ -93,7 +93,7 @@ def build_flowsheet():
m.fs.costing.add_annual_water_production(
m.fs.electroNP.properties_treated[0].flow_vol
)
m.fs.costing.add_LCOW(m.fs.electroNP.properties_treated[0].flow_vol)
m.fs.costing.add_LCOW(m.fs.AD.inlet.flow_vol[0])

# connections
m.fs.stream_adm1_translator = Arc(
Expand All @@ -105,7 +105,7 @@ def build_flowsheet():
pyo.TransformationFactory("network.expand_arcs").apply_to(m)

# Feed conditions based on mass balance in Flores-Alsina, where 0 terms are expressed as 1e-9
m.fs.AD.inlet.flow_vol.fix(
m.fs.AD.inlet.flow_vol[0].fix(
170 * units.m**3 / units.day
) # Double check this value
m.fs.AD.inlet.temperature.fix(308.15)
Expand Down Expand Up @@ -153,6 +153,9 @@ def build_flowsheet():
m.fs.electroNP.energy_electric_flow_mass.fix(0.044 * units.kWh / units.kg)
m.fs.electroNP.magnesium_chloride_dosage.fix(0.388)

# Costing
m.fs.costing.electroNP.phosphorus_recovery_value = 0

# scaling
for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
if "flow_vol" in var.name:
Expand All @@ -177,125 +180,101 @@ def build_flowsheet():
iscale.set_scaling_factor(var, 1e-1)
if "conc_mass_comp[X_I]" in var.name:
iscale.set_scaling_factor(var, 1e-1)
if "conc_mass_comp[S_O2]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[S_N2]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[S_NO3]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_H]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_PAO]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_AUT]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_MeOH]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_MeP]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_TSS]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[S_ch4]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_su]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_fa]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_c4]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_pro]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_ac]" in var.name:
iscale.set_scaling_factor(var, 1e4)
if "conc_mass_comp[X_h2]" in var.name:
iscale.set_scaling_factor(var, 1e4)

iscale.calculate_scaling_factors(m)

iscale.set_scaling_factor(m.fs.electroNP.properties_byproduct[0.0].flow_vol, 1e7)
iscale.set_scaling_factor(m.fs.AD.vapor_phase[0].pressure_sat, 1e-3)

skip_constraint = [
"S_O2",
"S_N2",
"X_AUT",
"X_PHA",
"S_I",
"S_NO3",
"S_Mg",
"X_MeP",
"X_MeOH",
"X_PAO",
"X_TSS",
"S_F",
"S_K",
"S_A",
"X_I",
"X_H",
"X_PP",
"X_S",
]

for i in skip_constraint:
m.fs.electroNP.solute_removal_equation[0.0, "Liq", i].deactivate()
m.fs.electroNP.solute_treated_equation[0.0, "Liq", i].deactivate()

m.fs.AD.initialize(outlvl=idaeslog.INFO_HIGH)
propagate_state(m.fs.stream_adm1_translator)
m.fs.translator_adm1_asm2d.initialize(outlvl=idaeslog.INFO_HIGH)
propagate_state(m.fs.stream_translator_electroNP)
m.fs.electroNP.initialize(outlvl=idaeslog.INFO_HIGH)
m.fs.costing.initialize()

solver = get_solver()

results = solver.solve(m, tee=True)
results = solve(m, tee=True)
return m, results


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_costing(m):
print("\nUnit Capital Costs\n")
print("\n----------Capital Cost----------")
total_capital_cost = value(
pyunits.convert(m.fs.costing.total_capital_cost, to_units=pyunits.USD_2018)
)
normalized_capex = total_capital_cost / value(
pyunits.convert(m.fs.AD.inlet.flow_vol[0], to_units=pyunits.m**3 / pyunits.hr)
)
print(f"Total Capital Costs: {total_capital_cost:.3f} $")
print(f"Normalized Capital Costs: {normalized_capex:.3f} $/m3/hr")
print("Capital Cost Breakdown")
for u in m.fs.costing._registered_unit_costing:
print(
u.name,
" : ",
value(units.convert(u.capital_cost, to_units=units.USD_2018)),
" : {price:0.3f} $".format(
price=value(pyunits.convert(u.capital_cost, to_units=pyunits.USD_2018))
),
)
print("\n----------Operation Cost----------")
total_operating_cost = value(
pyunits.convert(
m.fs.costing.total_operating_cost, to_units=pyunits.USD_2018 / pyunits.year
)
)
print(f"Total Operating Cost: {total_operating_cost:.3f} $/year")

opex_fraction = value(
pyunits.convert(
m.fs.costing.total_operating_cost, to_units=pyunits.USD_2018 / pyunits.year
)
/ pyunits.convert(
m.fs.AD.inlet.flow_vol[0], to_units=pyunits.m**3 / pyunits.year
)
/ m.fs.costing.LCOW
)
print(f"Operating cost fraction: {opex_fraction:.3f} $ opex / $ LCOW")

print("\nUtility Costs\n")
print("Operating Cost Breakdown")
for f in m.fs.costing.used_flows:
print(
f,
" : ",
value(
units.convert(
m.fs.costing.aggregate_flow_costs[f],
to_units=units.USD_2018 / units.year,
f.title(),
" : {price:0.3f} $/m3 feed".format(
price=value(
pyunits.convert(
m.fs.costing.aggregate_flow_costs[f],
to_units=pyunits.USD_2018 / pyunits.year,
)
/ pyunits.convert(
m.fs.AD.inlet.flow_vol[0],
to_units=pyunits.m**3 / pyunits.year,
)
)
),
)

print("")
total_capital_cost = value(
units.convert(m.fs.costing.total_capital_cost, to_units=units.USD_2018)
)
print(f"Total Capital Costs: {total_capital_cost:.2f} $")
total_operating_cost = value(
units.convert(
m.fs.costing.total_operating_cost, to_units=units.USD_2018 / units.year
)
)
print(f"Total Operating Costs: {total_operating_cost:.2f} $/year")
print("\n----------Energy----------")

electricity_intensity = value(
units.convert(
m.fs.electroNP.energy_electric_flow_mass, to_units=units.kWh / units.kg
pyunits.convert(
m.fs.costing.aggregate_flow_electricity / m.fs.AD.inlet.flow_vol[0],
to_units=units.kWh / units.m**3,
)
)
print(f"Electricity Intensity: {electricity_intensity:.4f} kWh/kg - P")
print(f"Electricity Intensity: {electricity_intensity:.4f} kWh/m3")

print("\n----------Levelized Cost----------")
LCOW = value(
units.convert(m.fs.costing.LCOW, to_units=units.USD_2018 / units.m**3)
pyunits.convert(m.fs.costing.LCOW, to_units=pyunits.USD_2018 / pyunits.m**3)
)
print(f"Levelized Cost of Water: {LCOW:.4f} $/m^3")
print(f"Levelized Cost of Water: {LCOW:.3f} $/m^3")


if __name__ == "__main__":
Expand All @@ -305,6 +284,7 @@ def display_costing(m):
{
"AD inlet": m.fs.AD.inlet,
"AD liquid outlet": m.fs.AD.liquid_outlet,
"AD vapor outlet": m.fs.AD.vapor_outlet,
"Translator outlet": m.fs.translator_adm1_asm2d.outlet,
"ElectroNP treated": m.fs.electroNP.treated,
"ElectroNP byproduct": m.fs.electroNP.byproduct,
Expand Down
Loading

0 comments on commit 06538b4

Please sign in to comment.