Simple dynamic solution running very slowly

[dhill2522]

I am having trouble formulating a dynamic model to run efficiently. I started with the lumped capacitance heat exchanger example and converted it to work for water in both streams. While the example problem runs quite quickly, my modified version will run for minutes on end with no solution.

Are there any well-known tips for efficient model formulation that I am missing? It seems that this should be a relatively simple problem, so I am not sure where things are going wrong.

My modified version of the problem is below:

import pyomo.environ as pe
from idaes.core import FlowsheetBlock
from idaes.models.unit_models import HeatExchangerLumpedCapacitance, HeatExchangerFlowPattern
from idaes.models.unit_models.heat_exchanger import delta_temperature_lmtd_callback
from idaes.models.properties import iapws95
from idaes.core.util.plot import plot_grid_dynamic
import matplotlib.pyplot as plt

m = pe.ConcreteModel()
m.fs = FlowsheetBlock(
    dynamic=True,
    time_set=[0, 300, 600, 900, 1200, 1500],
    time_units=pe.units.s
)

m.fs.prop_water = iapws95.Iapws95ParameterBlock()

m.fs.HE = HeatExchangerLumpedCapacitance(
    delta_temperature_callback=delta_temperature_lmtd_callback,
    cold_side_name="shell",
    hot_side_name="tube",
    shell={"property_package": m.fs.prop_water, "has_pressure_change": False},
    tube={"property_package": m.fs.prop_water, "has_pressure_change": False},
    flow_pattern=HeatExchangerFlowPattern.crossflow,
    dynamic_heat_balance=True,
    dynamic=False
)

m.discretizer = pe.TransformationFactory('dae.finite_difference')
m.discretizer.apply_to(m, nfe=10, wrt=m.fs.time, scheme="BACKWARD")

# Cold-side boundary conditions
shell_inlet_temperature = 288.15
shell_flow = 44004.14222
shell_area = 690073.9153
shell_hconv = 22
shell_pressure = 101325 # Pa
shell_enth = iapws95.htpx(T=shell_inlet_temperature*pe.units.K, P=shell_pressure*pe.units.Pa)
m.fs.HE.ua_cold_side[:].fix(shell_area * shell_hconv)
m.fs.HE.shell_inlet.flow_mol.fix(shell_flow)
m.fs.HE.shell_inlet.enth_mol[:].fix(shell_enth)
m.fs.HE.shell_inlet.pressure[:].fix(shell_pressure)

# Hot-side boundary conditions
tube_inlet_temperature = 384.35
tube_inlet_pressure = 7653000
tube_outlet_pressure = 7500000
tube_flow = 13896.84163
tube_area = 19542.2771
tube_hconv = 1000
m.fs.HE.ua_hot_side[:].fix(tube_area * tube_hconv)
m.fs.HE.tube_inlet.flow_mol[:].fix(tube_flow)
m.fs.HE.tube_inlet.pressure[:].fix(tube_inlet_pressure)

# number of tubes * tube mass * specific heat capacity
m.fs.HE.heat_capacity_wall = 1160 * 322 * 466
m.fs.HE.crossflow_factor.fix(0.8)

# Area has no effect on the result so long as it isn't zero
m.fs.HE.area.fix(1)

# Initialize the model with steady-state boundary conditions
m.fs.HE.initialize()
solver = pe.SolverFactory('ipopt')

# Activate the heat holdup term and add temperature disturbances
m.fs.HE.activate_dynamic_heat_eq()
for t in m.fs.time:
    enth = iapws95.htpx(T=288.15*pe.units.K, P=shell_pressure*pe.units.Pa)
    if 300 <= t < 600:
        enth = iapws95.htpx(T=(288.15 - 10)*pe.units.K, P=shell_pressure*pe.units.Pa)
    elif 600 <= t < 900:
        enth = iapws95.htpx(T=288.15*pe.units.K, P=shell_pressure*pe.units.Pa)
    elif 900 <= t < 1200:
        enth = iapws95.htpx(T=(288.15 + 10)*pe.units.K, P=shell_pressure*pe.units.Pa)
    elif t >= 1200:
        enth = iapws95.htpx(T=288.15*pe.units.K, P=shell_pressure*pe.units.Pa)
    m.fs.HE.shell_inlet.enth_mol[t].fix(enth)

print('starting solve')
solver.solve(m)

[andrewlee94]

@dhill2522 My first suggestion is to try using the new Diagnostics Toolbox to look for modeling issues:

Docs: https://idaes-pse.readthedocs.io/en/stable/reference_guides/core/util/diagnostics/diagnostics_toolbox.html#diagnostics-toolbox
Example: https://idaes-examples.readthedocs.io/en/latest/docs/diagnostics/diagnostics_toolbox_doc.html

There are many things that could be going wrong, and dynamic models are not easy to set up properly.

Two possibilities that immediately come to mind for me are:

1. Check that you have zero degrees of freedom. Whilst IDAES unit models support dynamic operation, you generally need to provide some additional constraints to define the outlet flowrate (e.g. flow driven or pressure driven operation).
2. The second thing to check is that the time steps you are providing are of a similar order of magnitude to the time constant of the system (which can be hard to know a priori). Having time steps that are too long (excessive error) or too short (stiffness) can result in a poorly posed model that is difficult to solve.

[dhill2522]

Thanks for the quick and very helpful response. This lead me to the correct solution, which was that I had neglected to specify the tube inlet enthalpy. Simple beginner mistake. I had not expected it to try to search for a solution if the problem was underspecified.