# Source code for idaes.examples.tutorials.Tutorial_3_Dynamic_Flowsheets

##############################################################################
# Institute for the Design of Advanced Energy Systems Process Systems
# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, by the
# software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia
#
# license information, respectively. Both files are also available online
# at the URL "https://github.com/IDAES/idaes-pse".
##############################################################################
"""
Demonstration and test flowsheet for a dynamic flowsheet.

"""

import matplotlib.pyplot as plt

# Import Pyomo libraries
from pyomo.environ import (ConcreteModel, SolverFactory, TransformationFactory,
Var, Constraint)
from pyomo.network import Arc

# Import IDAES core
from idaes.core import FlowsheetBlock

# Import Unit Model Modules
import idaes.property_models.examples.saponification_thermo as thermo_props
import idaes.property_models.examples.saponification_reactions as \
reaction_props

# Import Unit Model Modules
from idaes.unit_models import CSTR, Mixer

[docs]def main():
"""
Make the flowsheet object, fix some variables, and solve the problem
"""
# Create a Concrete Model as the top level object
m = ConcreteModel()

# Add a flowsheet object to the model
# time_set has points at 0 and 20 as the start and end of the domain,
# and a point at t=1 to allow for a step-change at this time
m.fs = FlowsheetBlock(default={"dynamic": True,
"time_set": [0, 1, 20]})

# Add property packages to flowsheet library
m.fs.thermo_params = thermo_props.SaponificationParameterBlock()
m.fs.reaction_params = reaction_props.SaponificationReactionParameterBlock(
default={"property_package": m.fs.thermo_params})

# Create unit models
m.fs.mix = Mixer(default={"dynamic": False,
"property_package": m.fs.thermo_params})
m.fs.Tank1 = CSTR(default={"property_package": m.fs.thermo_params,
"reaction_package": m.fs.reaction_params,
"has_holdup": True,
"has_equilibrium_reactions": False,
"has_heat_of_reaction": True,
"has_heat_transfer": True,
"has_pressure_change": False})
m.fs.Tank2 = CSTR(default={"property_package": m.fs.thermo_params,
"reaction_package": m.fs.reaction_params,
"has_holdup": True,
"has_equilibrium_reactions": False,
"has_heat_of_reaction": True,
"has_heat_transfer": True,
"has_pressure_change": False})

# Add pressure-flow constraints to Tank 1
m.fs.Tank1.height = Var(m.fs.time,
initialize=1.0,
doc="Depth of fluid in tank [m]")
m.fs.Tank1.area = Var(initialize=1.0,
doc="Cross-sectional area of tank [m^2]")
m.fs.Tank1.flow_coeff = Var(m.fs.time,
initialize=5e-5,
doc="Tank outlet flow coefficient")

def geometry(b, t):
return b.volume[t] == b.area*b.height[t]
m.fs.Tank1.geometry = Constraint(m.fs.time, rule=geometry)

def outlet_flowrate(b, t):
return b.control_volume.properties_out[t].flow_vol == \
b.flow_coeff[t]*b.height[t]**0.5
m.fs.Tank1.outlet_flowrate = Constraint(m.fs.time,
rule=outlet_flowrate)

# Add pressure-flow constraints to tank 2
m.fs.Tank2.height = Var(m.fs.time,
initialize=1.0,
doc="Depth of fluid in tank [m]")
m.fs.Tank2.area = Var(initialize=1.0,
doc="Cross-sectional area of tank [m^2]")
m.fs.Tank2.flow_coeff = Var(m.fs.time,
initialize=5e-5,
doc="Tank outlet flow coefficient")

m.fs.Tank2.geometry = Constraint(m.fs.time, rule=geometry)
m.fs.Tank2.outlet_flowrate = Constraint(m.fs.time,
rule=outlet_flowrate)

# Make Streams to connect units
m.fs.stream1 = Arc(source=m.fs.mix.outlet,
destination=m.fs.Tank1.inlet)

m.fs.stream2 = Arc(source=m.fs.Tank1.outlet,
destination=m.fs.Tank2.inlet)

# Discretize time domain
m.discretizer = TransformationFactory('dae.finite_difference')
m.discretizer.apply_to(m,
nfe=50,
wrt=m.fs.time,
scheme="BACKWARD")

TransformationFactory("network.expand_arcs").apply_to(m)

# Set inlet and operating conditions, and some initial conditions.
m.fs.mix.inlet_1.flow_vol.fix(0.5)
m.fs.mix.inlet_1.conc_mol_comp[:, "H2O"].fix(55388.0)
m.fs.mix.inlet_1.conc_mol_comp[:, "NaOH"].fix(100.0)
m.fs.mix.inlet_1.conc_mol_comp[:, "EthylAcetate"].fix(0.0)
m.fs.mix.inlet_1.conc_mol_comp[:, "SodiumAcetate"].fix(0.0)
m.fs.mix.inlet_1.conc_mol_comp[:, "Ethanol"].fix(0.0)
m.fs.mix.inlet_1.temperature.fix(303.15)
m.fs.mix.inlet_1.pressure.fix(101325.0)

m.fs.mix.inlet_2.flow_vol.fix(0.5)
m.fs.mix.inlet_2.conc_mol_comp[:, "H2O"].fix(55388.0)
m.fs.mix.inlet_2.conc_mol_comp[:, "NaOH"].fix(0.0)
m.fs.mix.inlet_2.conc_mol_comp[:, "EthylAcetate"].fix(100.0)
m.fs.mix.inlet_2.conc_mol_comp[:, "SodiumAcetate"].fix(0.0)
m.fs.mix.inlet_2.conc_mol_comp[:, "Ethanol"].fix(0.0)
m.fs.mix.inlet_2.temperature.fix(303.15)
m.fs.mix.inlet_2.pressure.fix(101325.0)

m.fs.Tank1.area.fix(1.0)
m.fs.Tank1.flow_coeff.fix(0.5)
m.fs.Tank1.heat_duty.fix(0.0)

m.fs.Tank2.area.fix(1.0)
m.fs.Tank2.flow_coeff.fix(0.5)
m.fs.Tank2.heat_duty.fix(0.0)

# Set initial conditions - accumulation = 0 at time = 0

# Initialize Units
m.fs.mix.initialize()

m.fs.Tank1.initialize(state_args={
"flow_vol": 1.0,
"conc_mol_comp": {"H2O": 55388.0,
"NaOH": 100.0,
"EthylAcetate": 100.0,
"SodiumAcetate": 0.0,
"Ethanol": 0.0},
"temperature": 303.15,
"pressure": 101325.0})

m.fs.Tank2.initialize(state_args={
"flow_vol": 1.0,
"conc_mol_comp": {"H2O": 55388.0,
"NaOH": 100.0,
"EthylAcetate": 100.0,
"SodiumAcetate": 0.0,
"Ethanol": 0.0},
"temperature": 303.15,
"pressure": 101325.0})

# Create a solver
solver = SolverFactory('ipopt')
results = solver.solve(m.fs)

# Make a step disturbance in feed and solve again
for t in m.fs.time:
if t >= 1.0:
m.fs.mix.inlet_2.conc_mol_comp[t, "EthylAcetate"].fix(90.0)
results = solver.solve(m.fs)

# Print results
print(results)

# For testing purposes
return(m, results)

def plot_results(m):
plt.figure("Tank 1 Outlet")
plt.plot(m.fs.time,
list(m.fs.Tank1.outlet.conc_mol_comp[:, "NaOH"].value),
label='NaOH')
plt.plot(m.fs.time,
list(m.fs.Tank1.outlet.conc_mol_comp[:, "EthylAcetate"].value),
label='EthylAcetate')
plt.plot(m.fs.time,
list(m.fs.Tank1.outlet.conc_mol_comp[:, "SodiumAcetate"].value),
label='SodiumAcetate')
plt.plot(m.fs.time,
list(m.fs.Tank1.outlet.conc_mol_comp[:, "Ethanol"].value),
label='Ethanol')
plt.legend()
plt.grid()
plt.xlabel("Time [s]")
plt.ylabel("Concentration [mol/m^3]")

plt.figure("Tank 2 Outlet")
plt.plot(m.fs.time,
list(m.fs.Tank2.outlet.conc_mol_comp[:, "NaOH"].value),
label='NaOH')
plt.plot(m.fs.time,
list(m.fs.Tank2.outlet.conc_mol_comp[:, "EthylAcetate"].value),
label='EthylAcetate')
plt.plot(m.fs.time,
list(m.fs.Tank2.outlet.conc_mol_comp[:, "SodiumAcetate"].value),
label='SodiumAcetate')
plt.plot(m.fs.time,
list(m.fs.Tank2.outlet.conc_mol_comp[:, "Ethanol"].value),
label='Ethanol')
plt.legend()
plt.grid()
plt.xlabel("Time [s]")
plt.ylabel("Concentration [mol/m^3]")
plt.show(block=True)

if __name__ == "__main__":
m, results = main()
plot_results(m)