Heat Exchanger 1D Unit Model with Two Property Packages¶
Problem Statement: In this example, we will be heating a benzene-toluene mixture using steam.
Tube Side Inlet
Flow Rate = 250 mol/s
Mole fraction (Benzene) = 0.4
Mole fraction (Toluene) = 0.6
Pressure = 101325 Pa
Temperature = 350 K
Shell Side Inlet
Flow Rate = 100 mol/s
Mole fraction (Steam) = 1
Pressure = 101325 Pa
Temperature = 450 K
This example will demonstrate the simulation of the 1D heat exchanger by fixing any 7 of the following degrees of freedom: - two of shell length, diameter, and area - two of tube length, diameter, and area - number of tubes - wall temperature (at all spatial points) - heat transfer coefficient (at all spatial points, for both shell and tube)
IDAES documentation reference for heat exchanger 1D model: https://idaes-pse.readthedocs.io/en/latest/reference_guides/model_libraries/generic/unit_models/heat_exchanger_1D.html
This example utilizes the simple 1D heat exchanger for a shell and tube
system. The IDAES library contains a more advanced ShellAndTube1D
heat exchanger supporting a 0D wall conduction model; more details on
the advanced 1D heat exchanger may be found
here.
Importing necessary tools¶
First, import the required IDAES and Pyomo modules. Note that the hotside (shell) and coldside (tube) properties leverage separate property packages:
# Import pyomo package
from pyomo.environ import ConcreteModel, Constraint, value, units
# Import idaes logger to set output levels
import idaes.logger as idaeslog
# Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model
from idaes.core import FlowsheetBlock
# Import the IAPWS property package to create a properties block for steam in the flowsheet
from idaes.models.properties import iapws95
from idaes.models.properties.iapws95 import htpx
from idaes.models.properties.modular_properties.base.generic_property import (
GenericParameterBlock)
# Import the BT_ideal property package to create a properties block for the tube side in the flowsheet
from idaes.models.properties.modular_properties.examples.BT_ideal \
import configuration
# Import the degrees_of_freedom function from the idaes.core.util.model_statistics package
from idaes.core.util.model_statistics import degrees_of_freedom
# Import the default IPOPT solver
from idaes.core.solvers import get_solver
# Import a heat exchanger unit
from idaes.models.unit_models.heat_exchanger_1D import (HeatExchanger1D,
HeatExchangerFlowPattern)
Setting up the flowsheet¶
Then, create the ConcreteModel foundation, attach the steady state
flowsheet, and declare the property parameter block that will used for
the shell and tube sides.
More information on this general workflow can be found here: https://idaes-pse.readthedocs.io/en/stable/how_to_guides/workflow/general.html
m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)
m.fs.properties_shell = iapws95.Iapws95ParameterBlock()
m.fs.properties_tube = GenericParameterBlock(**configuration)
In the 0D Heat Exchanger model example, geometry effects are ignored in favor of temperature gradients to calculate heat transfer. Here, we need to specify discretization, a flow configuration and wall conduction assumption. We specify the one-dimensional spatial discretization to use backwards finite difference approximations with 20 finite elements - these are the defaults if none are specified explicitly, and the shell and tube domains may be discretized differently if desired. The domains must use the same number of finite elements, as the exchanger is linear and the elements directly correspond according to the selected flow pattern.
The 1D Heat Exchanger supports the following flow configuration options:
HeatExchangerFlowPattern.cocurrent. Shell and tube flow in parallel, tube inlet transfers with shell inlet and similar with outlets, and temperature difference is greatest at the flow inlets (default).
HeatExchangerFlowPattern.countercurrent. Shell and tube flow in anti-parallel, tube inlet transfers with shell outlet and vice versa, and temperature difference changes minimally along the exchanger length.
The unit is created below:
# Create an instance of the heat exchanger unit, attaching it to the flowsheet
# Specify that the property packages to be used with the heat exchanger are the ones we created earlier.
m.fs.heat_exchanger = HeatExchanger1D(
hot_side_name="shell",
cold_side_name="tube",
shell={"property_package": m.fs.properties_shell,
"transformation_method": "dae.finite_difference",
"transformation_scheme": "BACKWARD"},
tube={"property_package": m.fs.properties_tube,
"transformation_method": "dae.finite_difference",
"transformation_scheme": "BACKWARD"},
finite_elements=20,
flow_type=HeatExchangerFlowPattern.cocurrent)
# Call the degrees_of_freedom function, get initial DOF
DOF_initial = degrees_of_freedom(m)
print("The initial DOF is {0}".format(DOF_initial))
The initial DOF is 31
assert DOF_initial == 31
Fixing input specifications¶
For this problem, we will fix the inlet conditions, re-evaluate the degrees of freedom to ensure the problem is square (i.e. DOF=0), and run two different options for unit specifications:
h = htpx(450*units.K, P = 101325*units.Pa) # calculate IAPWS enthalpy
m.fs.heat_exchanger.hot_side_inlet.flow_mol.fix(100) # mol/s
m.fs.heat_exchanger.hot_side_inlet.pressure.fix(101325) # Pa
m.fs.heat_exchanger.hot_side_inlet.enth_mol.fix(h) # J/mol
DOF_initial = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_initial))
The DOF is 28
m.fs.heat_exchanger.cold_side_inlet.flow_mol.fix(250) # mol/s
m.fs.heat_exchanger.cold_side_inlet.mole_frac_comp[0, "benzene"].fix(0.4)
m.fs.heat_exchanger.cold_side_inlet.mole_frac_comp[0, "toluene"].fix(0.6)
m.fs.heat_exchanger.cold_side_inlet.pressure.fix(101325) # Pa
m.fs.heat_exchanger.cold_side_inlet.temperature[0].fix(350) # K
DOF_final = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_final))
The DOF is 23
Option 1: Fix heat transfer coefficient (HTC) and dimensions of each domain¶
Below, we fix the heat exchanger area, length and heat transfer coefficient, which yields a fully defined problem for all finite elements with zero degrees of freedom that may be initialized and solved:
m.fs.heat_exchanger.area.fix(0.5) # m2
m.fs.heat_exchanger.length.fix(4.85) # m
m.fs.heat_exchanger.heat_transfer_coefficient.fix(500) # W/m2/K
DOF_final = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_final))
The DOF is 0
assert DOF_final == 0
Now that the problem is square (zero degress of freedom), we can initialize and solve the full model:
# Initialize the flowsheet, and set the output at INFO
m.fs.heat_exchanger.initialize(outlvl=idaeslog.INFO)
# Solve the simulation using ipopt
# Note: If the degrees of freedom = 0, we have a square problem
opt = get_solver()
solve_status = opt.solve(m, tee = True)
# Display a readable report
m.fs.heat_exchanger.report()
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.hot_side: Initialization Complete
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties: Starting initialization
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties: Dew and bubble point initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties: Equilibrium temperature initialization completed.
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties: State variable initialization completed.
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties: Phase equilibrium initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties: Property initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:01 [INFO] idaes.init.fs.heat_exchanger.cold_side: Initialization Complete
2023-03-04 01:48:02 [INFO] idaes.init.fs.heat_exchanger: Initialization Complete.
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200
**************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
**************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 2903
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 1121
Total number of variables............................: 907
variables with only lower bounds: 126
variables with lower and upper bounds: 308
variables with only upper bounds: 0
Total number of equality constraints.................: 907
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 7.51e+02 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 7.50e+00 8.47e-01 -1.0 4.26e-01 - 9.90e-01 9.90e-01h 1
2 0.0000000e+00 6.00e-02 2.00e+00 -1.0 3.47e-03 - 9.90e-01 9.92e-01h 1
3 0.0000000e+00 3.03e-08 9.98e+02 -1.0 2.77e-05 - 9.90e-01 1.00e+00h 1
Number of Iterations....: 3
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 3.0293449526652694e-08 3.0293449526652694e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 3.0293449526652694e-08 3.0293449526652694e-08
Number of objective function evaluations = 4
Number of objective gradient evaluations = 4
Number of equality constraint evaluations = 4
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 4
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 3
Total CPU secs in IPOPT (w/o function evaluations) = 0.049
Total CPU secs in NLP function evaluations = 0.095
EXIT: Optimal Solution Found.
====================================================================================
Unit : fs.heat_exchanger Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Area : 0.50000 : meter ** 2 : True : (None, None)
Length : 4.8500 : meter : True : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units shell Inlet shell Outlet tube Inlet tube Outlet
Molar Flow mole / second 100 100.00 - -
Mass Flow kilogram / second 1.8015 1.8015 - -
T kelvin 450.00 444.47 - -
P pascal 101325 1.0132e+05 - -
Vapor Fraction dimensionless 1.0000 1.0000 - -
Molar Enthalpy joule / mole 50977. 50780. - -
Total Molar Flowrate mole / second - - 250 250.00
Total Mole Fraction benzene dimensionless - - 0.40000 0.40000
Total Mole Fraction toluene dimensionless - - 0.60000 0.60000
Temperature kelvin - - 350 368.39
Pressure pascal - - 1.0132e+05 1.0132e+05
====================================================================================
from pyomo.environ import assert_optimal_termination
import pytest
# Check if termination condition is optimal
assert_optimal_termination(solve_status)
assert value(m.fs.heat_exchanger.hot_side_outlet.enth_mol[0]) == pytest.approx(
htpx(444.47*units.K, P = 101325*units.Pa), rel=1e-3)
assert value(m.fs.heat_exchanger.cold_side_outlet.temperature[0]) == pytest.approx(
368.39, rel=1e-3)
Option 2: Unfix shell length and fix shell outlet temperatures¶
In the previous case, we fixed the heat exchanger area, length and overall heat transfer coefficient. However, given that the models in IDAES are equation oriented, we can fix the outlet variables. For example, we can fix the hot outlet temperature and heat exchanger length, and solve for the heat exchanger area that will satisfy that condition.
m.fs.heat_exchanger.area.unfix()
m.fs.heat_exchanger.hot_side_outlet.enth_mol.fix(htpx(444.47*units.K, P = 101325*units.Pa))
DOF_final = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_final))
The DOF is 0
result = opt.solve(m)
print(result)
# Display a readable report
m.fs.heat_exchanger.report()
Problem:
- Lower bound: -inf
Upper bound: inf
Number of objectives: 1
Number of constraints: 907
Number of variables: 907
Sense: unknown
Solver:
- Status: ok
Message: Ipopt 3.13.2x3a Optimal Solution Found
Termination condition: optimal
Id: 0
Error rc: 0
Time: 0.15423965454101562
Solution:
- number of solutions: 0
number of solutions displayed: 0
====================================================================================
Unit : fs.heat_exchanger Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Area : 0.50043 : meter ** 2 : False : (None, None)
Length : 4.8500 : meter : True : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units shell Inlet shell Outlet tube Inlet tube Outlet
Molar Flow mole / second 100 100.00 - -
Mass Flow kilogram / second 1.8015 1.8015 - -
T kelvin 450.00 444.47 - -
P pascal 101325 1.0132e+05 - -
Vapor Fraction dimensionless 1.0000 1.0000 - -
Molar Enthalpy joule / mole 50977. 50780. - -
Total Molar Flowrate mole / second - - 250 250.00
Total Mole Fraction benzene dimensionless - - 0.40000 0.40000
Total Mole Fraction toluene dimensionless - - 0.60000 0.60000
Temperature kelvin - - 350 368.39
Pressure pascal - - 1.0132e+05 1.0132e+05
====================================================================================
# Check if termination condition is optimal
assert_optimal_termination(result)
assert value(m.fs.heat_exchanger.area) == pytest.approx(0.5, abs=1e-2)
assert value(m.fs.heat_exchanger.cold_side_outlet.temperature[0]) == pytest.approx(
368.39, rel=1e-3)