###############################################################################
# WaterTAP Copyright (c) 2021, 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/"
#
###############################################################################
from pyomo.environ import ConcreteModel, Objective, Expression, Constraint, TransformationFactory, value, Param, Block
from pyomo.environ import units as pyunits
from pyomo.network import Arc
from pyomo.util import infeasible
from idaes.core import FlowsheetBlock
from idaes.core.util.scaling import calculate_scaling_factors
from idaes.core.util.initialization import propagate_state
from watertap.examples.flowsheets.full_treatment_train.flowsheet_components import (pretreatment_NF,
desalination,
translator_block,
feed_block,
gypsum_saturation_index,
costing,
)
from watertap.examples.flowsheets.full_treatment_train.model_components import property_models
from watertap.examples.flowsheets.full_treatment_train.util import (solve_block,
check_dof,
check_build,
check_scaling)
"""Flowsheet example that satisfy minimum viable product requirements"""
[docs]def build_flowsheet_mvp_NF(m, **kwargs):
"""
Build a flowsheet with NF pretreatment and RO.
"""
# set up keyword arguments for the sections of treatment train
kwargs_pretreatment = {k: kwargs[k] for k in ('has_bypass', 'NF_type', 'NF_base',)}
kwargs_desalination = {k: kwargs[k] for k in ('has_desal_feed', 'is_twostage', 'has_ERD',
'RO_type', 'RO_base', 'RO_level',)}
# build flowsheet
property_models.build_prop(m, base='ion')
pretrt_port = pretreatment_NF.build_pretreatment_NF(m, **kwargs_pretreatment)
property_models.build_prop(m, base=kwargs['RO_base'])
desal_port = desalination.build_desalination(m, **kwargs_desalination)
property_models.build_prop(m, base='eNRTL')
translator_block.build_tb(m, base_inlet=kwargs['NF_base'],
base_outlet=kwargs['RO_base'],
name_str='tb_pretrt_to_desal')
# set up Arcs between pretreatment and desalination
m.fs.s_pretrt_tb = Arc(source=pretrt_port['out'], destination=m.fs.tb_pretrt_to_desal.inlet)
m.fs.s_tb_desal = Arc(source=m.fs.tb_pretrt_to_desal.outlet, destination=desal_port['in'])
# add gypsum saturation index calculations
gypsum_saturation_index.build(m, section='desalination', **kwargs_desalination)
gypsum_saturation_index.build(m, section='pretreatment', **kwargs_desalination)
# new initialization
if kwargs['NF_type'] == 'ZO':
m.fs.NF.area.fix(175)
if kwargs['has_bypass']:
m.fs.splitter.split_fraction[0, 'bypass'].fix(0.50)
m.fs.RO.area.fix(80)
m.fs.pump_RO.control_volume.properties_out[0].pressure.fix(60e5)
if kwargs['is_twostage']:
m.fs.RO2.area.fix(20)
m.fs.pump_RO2.control_volume.properties_out[0].pressure.fix(90e5)
# touch some properties used in optimization
if kwargs['is_twostage']:
product_water_sb = m.fs.mixer_permeate.mixed_state[0]
RO_waste_sb = m.fs.RO2.feed_side.properties[0, 1]
else:
product_water_sb = m.fs.RO.mixed_permeate[0]
RO_waste_sb = m.fs.RO.feed_side.properties[0, 1]
# NOTE: Building the costing here means it gets
# initialized during the simulation phase.
# This helps model stability.
m.fs.feed.properties[0].flow_vol
m.fs.feed.properties[0].conc_mol_phase_comp['Liq', 'Ca']
m.fs.tb_pretrt_to_desal.properties_in[0].flow_vol
m.fs.tb_pretrt_to_desal.properties_in[0].conc_mol_phase_comp['Liq', 'Ca']
product_water_sb.flow_vol
RO_waste_sb.flow_vol
m.fs.system_recovery = Expression(
expr=product_water_sb.flow_vol / m.fs.feed.properties[0].flow_vol)
m.fs.total_work = Expression(expr=m.fs.pump_RO.work_mechanical[0] +
(m.fs.pump_RO2.work_mechanical[0] if kwargs['is_twostage'] else 0.))
# annual water production
m.fs.treated_flow_vol = Expression(expr=product_water_sb.flow_vol)
costing.build_costing(m, **kwargs)
return m
def set_up_optimization(m, system_recovery=0.7, **kwargs_flowsheet):
is_twostage = kwargs_flowsheet['is_twostage']
# scale
calculate_scaling_factors(m)
# unfix variables
m.fs.splitter.split_fraction[0, 'bypass'].unfix()
m.fs.splitter.split_fraction[0, 'bypass'].setlb(0.001)
m.fs.splitter.split_fraction[0, 'bypass'].setub(0.99)
m.fs.NF.area.unfix()
m.fs.NF.area.setlb(10)
m.fs.NF.area.setub(1000)
m.fs.pump_RO.control_volume.properties_out[0].pressure.unfix()
m.fs.pump_RO.control_volume.properties_out[0].pressure.setlb(20e5)
m.fs.pump_RO.control_volume.properties_out[0].pressure.setub(75e5)
m.fs.RO.area.unfix()
m.fs.RO.area.setlb(10)
m.fs.RO.area.setub(300)
# Set lower bound for water flux at the RO outlet, based on a minimum net driving pressure, NDPmin
m.fs.RO.NDPmin = Param(initialize=1e5, mutable=True, units=pyunits.Pa)
m.fs.RO.flux_mass_phase_comp[0, 1, 'Liq', 'H2O'].setlb(value(m.fs.RO.A_comp[0, 'H2O']
* m.fs.RO.dens_solvent
* m.fs.RO.NDPmin))
if is_twostage:
m.fs.max_allowable_pressure = Param(initialize=120e5, mutable=True, units=pyunits.pascal)
m.fs.pump_RO2.control_volume.properties_out[0].pressure.unfix()
m.fs.pump_RO2.control_volume.properties_out[0].pressure.setlb(20e5)
m.fs.pump_RO2.control_volume.properties_out[0].pressure.setub(m.fs.max_allowable_pressure)
m.fs.RO2.area.unfix()
m.fs.RO2.area.setlb(10)
m.fs.RO2.area.setub(300)
# Set lower bound for water flux at the RO outlet, based on a minimum net driving pressure, NDPmin
m.fs.RO2.NDPmin = Param(initialize=1e5, mutable=True, units=pyunits.Pa)
m.fs.RO2.flux_mass_phase_comp[0, 1, 'Liq', 'H2O'].setlb(value(m.fs.RO2.A_comp[0, 'H2O']
* m.fs.RO2.dens_solvent
* m.fs.RO2.NDPmin))
# add additional constraints
# fixed system recovery
m.fs.system_recovery_target = Param(initialize=system_recovery, mutable=True)
m.fs.system_recovery_tol = Param(initialize=5e-3, mutable=True)
m.fs.eq_system_recovery = Constraint(
expr=(m.fs.system_recovery_target,
m.fs.system_recovery,
m.fs.system_recovery_target+m.fs.system_recovery_tol))
# saturation index
m.fs.max_saturation_index = Param(initialize=1.0, mutable=True)
m.fs.eq_max_saturation_index_desal = Constraint(
expr=m.fs.desal_saturation.saturation_index <= m.fs.max_saturation_index)
m.fs.max_conc_factor_target = Param(initialize=3.5, mutable=True)
m.fs.eq_max_conc_NF = Constraint(
expr=m.fs.NF.feed_side.properties_out[0].mass_frac_phase_comp['Liq', 'Ca']
<= m.fs.max_conc_factor_target * m.fs.feed.properties[0].mass_frac_phase_comp['Liq', 'Ca'])
# set objective
m.fs.objective = Objective(expr=m.fs.costing.LCOW)
# set additional constraints to limit local minima
# NOTE: doesn't seem necessary with new objective
if False:
m.fs.inequality_RO_area = Constraint(expr=m.fs.RO.area >= m.fs.RO2.area)
min_pressure_increase = 1e5
m.fs.inequality_RO_pressure = Constraint(
expr=m.fs.pump_RO.control_volume.properties_out[0].pressure + min_pressure_increase
<= m.fs.pump_RO2.control_volume.properties_out[0].pressure)
m.fs.inequality_RO_permeate = Constraint(
expr=m.fs.RO.permeate_side.properties_mixed[0].flow_vol_phase['Liq']
>= m.fs.RO2.permeate_side.properties_mixed[0].flow_vol_phase['Liq'])
check_dof(m, dof_expected=6 if is_twostage else 4)
# solve_block(m, tee=False, fail_flag=True)
def optimize(m):
solve_block(m, tee=False, fail_flag=True)
def solve_flowsheet_mvp_NF(**kwargs):
m = ConcreteModel()
m.fs = FlowsheetBlock(default={"dynamic": False})
build_flowsheet_mvp_NF(m, **kwargs)
TransformationFactory("network.expand_arcs").apply_to(m)
# scale
pretreatment_NF.scale_pretreatment_NF(m, **kwargs)
calculate_scaling_factors(m.fs.tb_pretrt_to_desal)
desalination.scale_desalination(m, **kwargs)
calculate_scaling_factors(m)
# initialize
optarg = {'nlp_scaling_method': 'user-scaling'}
pretreatment_NF.initialize_pretreatment_NF(m, **kwargs)
m.fs.pretrt_saturation.properties.initialize(optarg=optarg)
propagate_state(m.fs.s_pretrt_tb)
m.fs.tb_pretrt_to_desal.initialize(optarg=optarg)
propagate_state(m.fs.s_tb_desal)
desalination.initialize_desalination(m, **kwargs)
m.fs.desal_saturation.properties.initialize()
m.fs.costing.initialize()
# check_build(m)
# check_scaling(m)
check_dof(m)
solve_block(m, tee=False, fail_flag=True)
pretreatment_NF.display_pretreatment_NF(m, **kwargs)
m.fs.tb_pretrt_to_desal.report()
desalination.display_desalination(m, **kwargs)
print('desalination solubility index:', value(m.fs.desal_saturation.saturation_index))
print('pretreatment solubility index:', value(m.fs.pretrt_saturation.saturation_index))
print('water recovery:', value(m.fs.system_recovery))
print('LCOW:', value(m.fs.costing.LCOW))
print('CP modulus:', value(m.fs.desal_saturation.cp_modulus))
return m
def solve_optimization(system_recovery=0.75, **kwargs_flowsheet):
m = solve_flowsheet_mvp_NF(**kwargs_flowsheet)
print('\n****** Optimization *****\n')
set_up_optimization(m, system_recovery=system_recovery, **kwargs_flowsheet)
optimize(m)
pretreatment_NF.display_pretreatment_NF(m, **kwargs_flowsheet)
m.fs.tb_pretrt_to_desal.report()
desalination.display_desalination(m, **kwargs_flowsheet)
costing.display_costing(m)
print('desalination saturation index:', value(m.fs.desal_saturation.saturation_index))
print('pretreatment saturation index:', value(m.fs.pretrt_saturation.saturation_index))
print('pretreatment Ca concentration factor:', value(
m.fs.NF.feed_side.properties_out[0].mass_frac_phase_comp['Liq', 'Ca'] /
m.fs.feed.properties[0].mass_frac_phase_comp['Liq', 'Ca']))
print('water recovery:', value(m.fs.system_recovery))
print('CP modulus:', value(m.fs.desal_saturation.cp_modulus))
return m
if __name__ == "__main__":
import sys
kwargs_flowsheet = {
'has_bypass': True, 'has_desal_feed': False, 'is_twostage': True, 'has_ERD': True,
'NF_type': 'ZO', 'NF_base': 'ion',
'RO_type': '0D', 'RO_base': 'TDS', 'RO_level': 'detailed'}
if len(sys.argv) == 1:
m = solve_flowsheet_mvp_NF(**kwargs_flowsheet)
else:
m = solve_optimization(system_recovery=float(sys.argv[1]), **kwargs_flowsheet)