Initialization Failure in Heat Exchanger for Waste Heat Recovery

[Hanli131]

Description
To simplify the problem, I built a very basic system as shown below: The outlet gas from a 0D heat exchanger is heated by a heater (representing a Solid Oxide Fuel Cell in the actual system) and then used to preheat the inlet gas. However, the initialization process fails, and I don't know what the issue is.

The system fails during initialization when the heater outlet is connected to the heat exchanger hot-side inlet.
If I remove this connection (m.fs.s03), the model initializes and runs successfully.
The model degrees of freedom (DOF) seem reasonable after fixing necessary variables.

Questions
What could be causing the initialization failure in this setup?
Which resources or tutorials should I refer to for troubleshooting such issues?

Any insights or suggestions would be greatly appreciated. Thanks in advance!

Here is the code:

# Import required packages
from idaes.models_extra.power_generation.properties.natural_gas_PR import get_prop
from pyomo.environ import ConcreteModel, SolverFactory, value, units, TransformationFactory, Expression
from idaes.core import FlowsheetBlock
import idaes.logger as idaeslog
from idaes.models.properties import iapws95
from idaes.models.unit_models import Feed
from idaes.models.unit_models.heater import Heater
from idaes.models.unit_models.heat_exchanger import HeatExchanger, delta_temperature_amtd_callback
from pyomo.network import Arc, SequentialDecomposition
from idaes.core.util.model_statistics import degrees_of_freedom

# Create the ConcreteModel and the FlowsheetBlock, and attach the flowsheet block to it.
m = ConcreteModel()

m.fs = FlowsheetBlock(dynamic=False)

m.fs.properties = iapws95.Iapws95ParameterBlock()

m.fs.feed = Feed(property_package=m.fs.properties)
m.fs.heater = Heater(property_package=m.fs.properties)
m.fs.heat_exchanger = HeatExchanger(
        delta_temperature_callback=delta_temperature_amtd_callback,
        hot_side_name="shell",
        cold_side_name="tube",
        shell={"property_package": m.fs.properties},
        tube={"property_package": m.fs.properties})

m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.heat_exchanger.cold_side_inlet)
m.fs.s02 = Arc(source=m.fs.heat_exchanger.cold_side_outlet, destination=m.fs.heater.inlet)
m.fs.s03 = Arc(source=m.fs.heater.outlet, destination=m.fs.heat_exchanger.hot_side_inlet)

TransformationFactory("network.expand_arcs").apply_to(m)
DOF_initial = degrees_of_freedom(m)
print("The initial DOF is {0}".format(DOF_initial))

# Fix the stream inlet conditions
m.fs.feed.flow_mol[0].fix(100) #mol/s
m.fs.feed.pressure[0].fix(101325) #Pa
m.fs.feed.enth_mol[0].fix(value(iapws95.htpx(T=293*units.K, P=101325*units.Pa)))


m.fs.heat_exchanger.overall_heat_transfer_coefficient[0].fix(500) # W/m2/K
m.fs.heat_exchanger.area.fix(50)
#m.fs.heat_exchanger.shell_inlet.flow_mol.fix(100) # mol/s
#m.fs.heat_exchanger.shell_inlet.enth_mol[0].fix(value(iapws95.htpx(T=1073*units.K, P=101325*1.5*units.Pa)))
#m.fs.heat_exchanger.shell_inlet.pressure[0].fix(101325) #Pa

m.fs.heater.outlet.enth_mol.fix(value(iapws95.htpx(T=1073*units.K, P=101325*units.Pa)))

DOF_initial = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_initial))

# Initialize the flowsheet, and set the logger level to INFO
seq = SequentialDecomposition()
seq.options.select_tear_method = "heuristic"
seq.options.tear_method = "Wegstein"
seq.options.iterLim = 5

def function(unit):
        unit.initialize(outlvl=idaeslog.INFO)
        
seq.run(m, function)

# Solve the model
from idaes.core.solvers import get_solver
solver = get_solver()


results = solver.solve(m, tee=True)
m.fs.heat_exchanger.report()

[dallan-keylogic]

Running m.fs.heat_exchanger.report() reveals a problem rather quickly:

====================================================================================
Unit : fs.heat_exchanger                                                   Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables:

    Key            : Value   : Units                           : Fixed : Bounds
        Delta T In :  788.33 :                          kelvin : False : (None, None)
       Delta T Out : -782.33 :                          kelvin : False : (None, None)
           HX Area :  50.000 :                      meter ** 2 :  True : (0, None)
    HX Coefficient :  500.00 : kilogram / kelvin / second ** 3 :  True : (0, None)
         Heat Duty :  74941. :                            watt : False : (None, None)

    Expressions:

    Key             : Value  : Units
    Delta T Driving : 2.9976 : kelvin

------------------------------------------------------------------------------------
    Stream Table
                         Units         shell Inlet  shell Outlet  tube Inlet  tube Outlet
    Molar Flow          mole / second      1.0000        1.0000       100.00      100.00
    Mass Flow       kilogram / second    0.018015      0.018015       1.8015      1.8015
    T                          kelvin      1093.4        270.49       293.00      302.95
    P                          pascal  1.1032e+07    1.1032e+07   1.0132e+05  1.0132e+05
    Vapor Fraction      dimensionless      1.0000        0.0000       0.0000      0.0000
    Molar Enthalpy       joule / mole      74941.      0.011021       1502.1      2251.5
====================================================================================

The shell flow rate is miniscule (1 mol/s) compared to the tube flow rate (100 mol/s). The shell side temperature crashes into the property package's lower bound for temperature, causing the infeasibility of the initialization.

So why is the shell flow rate so small? Let's look at the output from the sequential decomposition:

2025-03-14 16:37:46 [INFO] idaes.init.fs.feed: Initialization Complete.
2025-03-14 16:37:46 [INFO] idaes.init.fs.heater.control_volume: Initialization Complete
2025-03-14 16:37:46 [INFO] idaes.init.fs.heater: Initialization Complete: optimal - Optimal Solution Found
2025-03-14 16:37:46 [INFO] idaes.init.fs.heat_exchanger.hot_side: Initialization Complete
2025-03-14 16:37:46 [INFO] idaes.init.fs.heat_exchanger.cold_side: Initialization Complete
2025-03-14 16:37:47 [INFO] idaes.init.fs.heat_exchanger: Initialization Completed, infeasible - Converged to a locally infeasible point. Problem may be infeasible.

It is initializing the heater before the heat exchanger---the tear stream must be s2. However, it has no initial guess for any of the state variables in the heater's inlet, so it uses the values at which those variables were initialized---in this case 1 mol/s for flow rate.

The solution is then to provide a better initial guess at the tear stream. It appears Pyomo's SequentialDecomposition supports providing initial guesses, but I haven't used it in that fashion. Instead, you could create a tear stream manually by deactivating an arc constraint, provide an initial condition, solve the flowsheet, then reactivate the arc constraint for a final solve:

# Import required packages
from idaes.models_extra.power_generation.properties.natural_gas_PR import get_prop
from pyomo.environ import ConcreteModel, SolverFactory, value, units, TransformationFactory, Expression
from idaes.core import FlowsheetBlock
from idaes.core.util.initialization import propagate_state
import idaes.logger as idaeslog
from idaes.models.properties import iapws95
from idaes.models.unit_models import Feed
from idaes.models.unit_models.heater import Heater
from idaes.models.unit_models.heat_exchanger import HeatExchanger, delta_temperature_amtd_callback
from pyomo.network import Arc, SequentialDecomposition
from idaes.core.util.model_statistics import degrees_of_freedom

# Create the ConcreteModel and the FlowsheetBlock, and attach the flowsheet block to it.
m = ConcreteModel()

m.fs = FlowsheetBlock(dynamic=False)

m.fs.properties = iapws95.Iapws95ParameterBlock()

m.fs.feed = Feed(property_package=m.fs.properties)
m.fs.heater = Heater(property_package=m.fs.properties)
m.fs.heat_exchanger = HeatExchanger(
        delta_temperature_callback=delta_temperature_amtd_callback,
        hot_side_name="shell",
        cold_side_name="tube",
        shell={"property_package": m.fs.properties},
        tube={"property_package": m.fs.properties})

m.fs.s01 = Arc(source=m.fs.feed.outlet, destination=m.fs.heat_exchanger.cold_side_inlet)
m.fs.s02 = Arc(source=m.fs.heat_exchanger.cold_side_outlet, destination=m.fs.heater.inlet)
m.fs.s03 = Arc(source=m.fs.heater.outlet, destination=m.fs.heat_exchanger.hot_side_inlet)

TransformationFactory("network.expand_arcs").apply_to(m)
DOF_initial = degrees_of_freedom(m)
print("The initial DOF is {0}".format(DOF_initial))

# Fix the stream inlet conditions
m.fs.feed.flow_mol[0].fix(100) #mol/s
m.fs.feed.pressure[0].fix(101325) #Pa
m.fs.feed.enth_mol[0].fix(value(iapws95.htpx(T=293*units.K, P=101325*units.Pa)))


m.fs.heat_exchanger.overall_heat_transfer_coefficient[0].fix(500) # W/m2/K
m.fs.heat_exchanger.area.fix(50)


m.fs.heater.outlet.enth_mol.fix(value(iapws95.htpx(T=1073*units.K, P=101325*units.Pa)))

DOF_initial = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_initial))

# Provide initial guess for the shell inlet and create tear stream
m.fs.heat_exchanger.shell_inlet.flow_mol.fix(100) # mol/s
m.fs.heat_exchanger.shell_inlet.enth_mol[0].fix(value(iapws95.htpx(T=1073*units.K, P=101325*1.5*units.Pa)))
m.fs.heat_exchanger.shell_inlet.pressure[0].fix(101325) #Pa
m.fs.s03_expanded.deactivate()

DOF_initial = degrees_of_freedom(m)
print("The DOF is {0} after creating tear stream".format(DOF_initial))

m.fs.feed.initialize(outlvl=idaeslog.INFO)
propagate_state(m.fs.s01)
m.fs.heat_exchanger.initialize(outlvl=idaeslog.INFO)    
propagate_state(m.fs.s02)
m.fs.heater.initialize(outlvl=idaeslog.INFO)

# Solve the model
from idaes.core.solvers import get_solver
solver = get_solver()
results = solver.solve(m, tee=True)

# Reactivate tear stream, and unfix shell side initial conditions
m.fs.heat_exchanger.shell_inlet.flow_mol.unfix()
m.fs.heat_exchanger.shell_inlet.enth_mol[0].unfix()
m.fs.heat_exchanger.shell_inlet.pressure[0].unfix()
m.fs.s03_expanded.activate()

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

m.fs.heat_exchanger.report()

For more advanced diagnostics, check out the Diagnostics Toolbox.

[Hanli131]

Thank you very much for your prompt response! I’ll review the materials and make corrections soon! ：）

[Hanli131]

Sorry to bother you again. Could you provide an example of solving this issue by setting an initial guess at the tear stream?

I tried adding the following code to the source, and the output confirmed that fs.s02 is the tear stream. However, the initialization still failed.

When I ran m.fs.heat_exchanger.report(), I found that the initial guess did not take effect—the hot-side flow rate of the heat exchanger remained at 1 mol/s.

for o in heuristic_tear_set:
    print(o.name)

seq.set_guesses_for(m.fs.heater.inlet, {
    m.fs.heater.inlet.flow_mol: 100,  # Initial flow guess in mol/s
})

[dallan-keylogic]

It looks like the issue is that you're treating m.fs.heater.inlet.flow_mol like an unindexed component, when it is in fact indexed by time.

The guesses should be a dict that maps the following:

Port Member Name -> Value

Or, for indexed members, multiple dicts that map:

Port Member Name -> Index -> Value

I would guess that this looks like

seq.set_guesses_for(
    m.fs.heater.inlet,
    {
        m.fs.heater.inlet.flow_mol: {
            0: 100,    # Initial flow guess in mol/s
        }
    }
)

If the SequentialDecomposition didn't give some sort of warning for your input like Ignoring input xyz, you should open an issue on Pyomo about it.

[Hanli131]

Thank you for your response. Referring to your suggestion, after trying various initialization approaches, I finally managed to get it working with the following code.

seq.set_guesses_for(m.fs.heater.inlet, {
                    "flow_mol":{0:100}  # Initial flow guess in mol/s
                    })

It seems that using the generic name “flow_mol” (instead of “m.fs.heater.inlet.flow_mol”) is sufficient for stream guessing.

Once again, I really appreciate your help! 💐💐